Logical Methods in Computer Science 
Vol. 4 (4:6) 2008, pp. 1-43 
www.lmcs-online.org 



Submitted Aug. 31 , 2007 
Published Nov. 11,2008 



FLOW FASTER: EFFICIENT DECISION ALGORITHMS FOR 
PROBABILISTIC SIMULATIONS * 

LIJUN ZHANG a , HOLGER HERMANNS 6 , FRIEDRICH EISENBRAND c , AND DAVID N. JANSEN d 

a ' b Department of Computer Science, Saarland University, Germany 
e-mail address: {zhang,hermanns}@cs.uni-sb.de 

c Department of Mathematics, EPFL, Switzerland 
e-mail address: friedrich.eisenbrand@epfl.ch 

d Software Modelling and Verification, RWTH Aachen University, Germany, and Model-Based Sys- 
tem Development, Radboud University, Nijmegen, The Netherlands 
e-mail address: D.Jansen@cs.ru.nl 



Abstract. Strong and weak simulation relations have been proposed for Markov chains, 
while strong simulation and strong probabilistic simulation relations have been proposed 
for probabilistic automata. This paper investigates whether they can be used as effec- 
tively as their non-probabilistic counterparts. It presents drastically improved algorithms 
to decide whether some (discrete- or continuous-time) Markov chain strongly or weakly 
simulates another, or whether a probabilistic automaton strongly simulates another. The 
key innovation is the use of parametric maximum flow techniques to amortize compu- 
tations. We also present a novel algorithm for deciding strong probabilistic simulation 
preorders on probabilistic automata, which has polynomial complexity via a reduction 
to an LP problem. When extending the algorithms for probabilistic automata to their 
continuous-time counterpart, we retain the same complexity for both strong and strong 
probabilistic simulations. 



Many verification methods have been introduced to prove the correctness of systems ex- 
ploiting rigorous mathematical foundations. As one of the automatic verification techniques, 
model checking has successfully been applied to automatically find errors in complex sys- 
tems. The power of model checking is limited by the state space explosion problem. Notably, 
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minimizing the system to the bisimulation |34[ [35] quotient is a favorable approach. As a 
more aggressive attack to the problem, simulation relations |33| have been proposed for 
these models. In particular, they provide the principal ingredients to perform abstrac- 
tions of the models, while preserving safe CTL properties (formulas with positive universal 
path-quantifiers only) [16j . 

Simulation relations are preorders on the state space such that whenever state s' sim- 
ulates state s (written s ^ s') then s' can mimic all stepwise behaviour of s, but s' may 
perform steps that cannot be matched by s. One of the interesting aspects of simulation 
relations is that they allow a verification by "local" reasoning. Based on this, efficient 
algorithms for deciding simulation preorders have been proposed in [10\ 123]. 

Randomisation has been employed widely for performance and dependability models, 
and consequently the study of verification techniques of probabilistic systems with and with- 
out nondeterminism has drawn a lot of attention in recent years. A variety of equivalence 
and preorder relations, including strong and weak simulation relations, have been introduced 
and widely considered for probabilistic models. In this paper we consider discrete-time 
Markov chains (DTMCs) and discrete-time probabilistic automata (PAs) [39]. PAs extend 
labelled transition systems (LTSs) with probabilistic selection, or, viewed differently, extend 
DTMCs with nondeterminism. They constitute a natural model of concurrent computation 
involving random phenomena. In a PA, a labelled transition leads to a probability distri- 
bution over the set of states, rather than a single state. The resulting model thus exhibits 
both non-deterministic choice (as in LTSs) and probabilistic choice (as in Markov chains). 

Strong simulation relations have been introduced [26} [30] for probabilistic systems. 
For s ^ s' (s' strongly simulates s), it is required that every successor distribution of s 
via action a (called a-successor distribution) has a corresponding a-successor distribution 
at s' . Correspondence of distributions is naturally defined with the concept of weight 
functions [26]. In the context of model checking, strong simulation relations preserve safe 
PCTL formulas [40]. Probabilistic simulation [40] is a relaxation of strong simulation in the 
sense that it allows for convex combinations of multiple distributions belonging to equally 
labelled transitions. More concretely, it may happen that for an a-successor distribution fi 
of s, there is no a-successor distribution of s' which can be related to fj,, yet there exists a 
so-called a- combined transition, a convex combination of several a-successor distributions 
of s'. Probabilistic simulation accounts for this and is thus coarser than strong simulation, 
but still preserves the same class of PCTL-properties as strong simulation does. 

Apart from discrete time models, this paper considers continuous-time Markov chains 
(CTMCs) and continuous-time probabilistic automata (CPAs) [29\ 142] . In CPAs, the tran- 
sition delays are governed by exponential distributions. CPAs can be considered also as 
extensions of CTMCs with nondeterminism. CPAs are natural foundational models for 
various performance and dependability modelling formalisms including stochastic activity 
networks |37j . generalised stochastic Petri nets ]32] and interactive Markov chains |24j . 
Strong simulation and probabilistic simulation have been introduced for continuous-time 
models [91H3]. For CTMCs, s ^ s' requires that s ^ s' holds in the embedded DTMC, and 
additionally, state s' must be "faster" than s which manifests itself by a higher exit rate. 
Both strong simulation and probabilistic simulation preserve safe CSL formulas [3], which 
is a continuous stochastic extension of PCTL, tailored to continuous-time models. 

Weak simulation is proposed in [9] for Markov chains. In weak simulation, the successor 
states are split into visible and invisible parts, and the weight function conditions are 
only imposed on the transitions leading to the visible parts of the successor states. Weak 
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simulation is strictly coarser than the afore- mentioned strong simulation for Markov chains, 
thus allows further reduction of the state space. It preserves the safe PCTL- and CSL- 
properties without the next state formulas for DTMCs and CTMCs respectively [9]. 



Decision algorithms for strong and weak simulations over Markov chains, and for strong 
simulation over probabilistic automata are not efficient, which makes it as yet unclear 
whether they can be used as effectively as their non-probabilistic counterparts. In this paper 
we improve efficient decision algorithms, and devise new algorithms for deciding strong and 
strong probabilistic simulations for probabilistic automata. Given the simulation preorder, 
the simulation quotient automaton is in general smaller than the bisimulation quotient 
automaton. Then, for safety and liveness properties, model checking can be performed 
on this smaller quotient automata. The study of decision algorithms is also important 
for specification relations: The model satisfies the specification if the automaton for the 
specification simulates the automaton for the model. In many applications the specification 
cannot be easily expressed by logical formulas: it is rather a probabilistic model itself. 
Examples of this kind include various recent wireless network protocols, such as ZigBee |21| . 
Firewire Zeroconf [11] . or the novel IEEE 802. lie, where the central mechanism is selecting 
among different-sided dies, readily expressible as a probabilistic automaton [3T] . 

The common strategy used by decision algorithms for simulations is as follows. The 
algorithm starts with a relation R which is guaranteed to be coarser than the simulation pre- 
order ^. Then, the relation R is successively be refined. In each iteration of the refinement 
loop, pairs (s, s') are eliminated from the relation R if the corresponding simulation condi- 
tions are violated with respect to the current relation. In the context of labelled transitions 
systems, this happens if s has a successor state t, but we cannot find a successor state t' of s' 
such that (t, t') is also in the current relation R. For DTMCs, this correspondence is formu- 
lated by the existence of a weight function for distributions (P(s, -),P(s', •)) with respect 
to the current relation R. Checking this weight function condition amounts to checking 
whether there is a maximum flow over the network constructed out of (P(s, •), P(V, •)) and 
the current relation R. The complexity for one such check is however rather expensive, it 
has time complexity C(n 3 /logn). If the iterative algorithm reaches a fix-point, the strong 
simulation preorder is obtained. The number of iterations of the refinement loop is at most 
0(n 2 ), and the overall complexity [3] amounts to C(n 7 /logn) in time and 0(n 2 ) in space. 

Fixing a pair (s, s'), we observe that the networks for this pair across iterations of the 
refinement loop are very similar: They differ from iteration to iteration only by deletion of 
some edges induced by the successive clean up of R. We exploit this by adapting a para- 
metric maximum flow algorithm [18J to solve the maximum flow problems for the arising 
sequences of similar networks, hence arriving at efficient simulation decision algorithms. 
The basic idea is that all computations concerning the pair (s, s') can be performed in an 
incremental way: after each iteration we save the current network together with maximum 
flow information. Then, in the next iteration, we update the network, and derive the max- 
imum flow while using the previous maximum flow function. The maximum flow problems 
for the arising sequences of similar networks with respect to the pair (s, s') can be computed 
in time C(|F| 3 ) where \V\ is the number of nodes of the network. This leads to an overall 
time complexity 0(m 2 n) for deciding the simulation preorder. Because of the storage of 
the networks, the space complexity is increased to C(m 2 ). Especially in the very common 
case where the state fanout of a model is bounded by a constant g (and hence m < gn), 
our strong simulation algorithm has time and space complexity 0(n 2 ). The algorithm can 
be extended easily to handle CTMCs with same time and space complexity. For weak 
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simulation on Markov chains, the parametric maximum flow technique cannot be applied 
directly. Nevertheless, we manage to incorporate the parametric maximum flow idea into a 
decision algorithm with time complexity C(m 2 n 3 ) and space complexity 0(n 2 ). An earlier 
algorithm [6] uses LP problems |27l [38] as subroutines. The maximum flow problem is a 
special instance of an LP problem but can be solved much more efficiently pQ. 

We extend the algorithm to compute strong simulation preorder to also work on PAs. 
It takes the skeleton of the algorithm for Markov chains: It starts with a relation R which is 
coarser than ^, and then refines R until ^ is achieved. In the refinement loop, a pair (s, s') 
is eliminated if the corresponding simulation conditions are violated with respect to the 
current relation. For PAs, this means that there exists an a-successor distribution /i of s, 
such that for all a-successor distribution /i' of s', we cannot find a weight function for (/j,, //) 
with respect to the current relation R. Again, as for Markov chains, the existence of such 
weight functions can be reduced to maximum flow problems. Combining with the parametric 
maximum flow algorithm |18| . we arrive at the same time complexity 0(m 2 n) and space 
complexity 0{m?) as for Markov chains. The above maximum flow based procedure cannot 
be applied to deal with strong probabilistic simulation for PAs. The reason is that an a- 
combined transition of state s is a convex combination of several a-successor distributions of 
s, thus induces uncountable many such possible combined transitions. The computational 
complexity of deciding strong probabilistic simulation has not been investigated before. We 
show that it can be reduced to solving LP problems. The idea is that we introduce for each 
a-successor distribution a variable, and then reformulates the requirements concerning the 
combined transitions by linear constraints over these variables. This allows us to construct 
a set of LP problem such that whether a pair (s, s') should be thrown out of the current 
pair R is equivalent to whether each of the LP problem has a solution. 

The algorithms for PAs are then extended to handle their continuous-time analogon, 
CPAs. In the algorithm, for each pair (s, s') in the refinement loop, an additional rate 
condition is ensured by an additional check via comparing the appropriate rates of s and 
s' . The resulting algorithm has the same time and space complexity. 

Related Works. In the non-probabilistic setting, the most efficient algorithms for deciding 
simulation preorders have been proposed in [10, 23J. The complexity is 0(mn) where n 
and m denote the number of states and transitions of the transition system respectively. 
For Markov chains, Derisavi et al. [T7] presented an O(mlogn) algorithm for strong bisim- 
ulation. Weak bisimulation for DTMCs can be computed in C(n 3 ) time [5J. For strong 
simulation, Baier et al. [3J introduced a polynomial decision algorithm with complexity 
C(n 7 /logn), by tailoring a network flow algorithm [20] to the problem, embedded into an 
iterative refinement loop. In [6j, Baier et al. proved that weak simulation is decidable in 
polynomial time by reducing it to linear programming (LP) problems. For a subclass of PAs 
(reactive systems), Huynh and Tian [25] presented an O(mlogn) algorithm for computing 
strong bisimulation. Cattani and Segala [12] have presented decision algorithms for strong 
and ^simulation for PAs. They reduced the decision problems to LP problems. To compute 
the coarsest strong simulation for PAs, Baier et al. [3] presented an algorithm which reduces 
the query whether a state strongly simulates another to a maximum flow problem. Their 
algorithm has complexity 0((mn 6 + m 2 n 3 )/ log njH Recently, algorithm for computing 
simulation and bisimulation metrics for concurrent games [13j has been studied. 

^The m used in paper !3 S is slightly different from the m as we use it. A detailed comparison is provided 
later, in Remark 14. Ill of Section f4. 31 
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Outline of The Paper. The paper proceeds by recalling the definition of the models and 
simulation relations in Section [2j In Section [3] we give a short interlude on maximum flow 
problems. In Section 2] we present a combinatorial method to decide strong simulations. 
In this section we also introduce new decision algorithms for deciding strong probabilistic 
simulations for PAs and CPAs. In Section [5] we focus on algorithms for weak simulations. 
Section [6] concludes the paper. 

2. Preliminaries 

In Subsection 12.11 we recall the definitions of fully probabilistic systems, discrete- and 
continuous-time Markov chains [41], and the nondeterministic extensions of these discrete- 
time [40J and continuous-time models [36l E] . In Subsection 12.21 we recall the definition of 
simulation relations. 

2.1. Markov Models. Firstly, we introduce some general notations. Let X, Y be finite 
sets. For / : X -> R, let f(A) denote Y, X &A for all AC X. For / : X x Y -> R, we let 
f(x, A) denote Yly^A /( x > V) f° r ah ^ £ A and A CY, and f(A, y) is defined similarly. Let 
AP be a fixed, finite set of atomic propositions. 

For a finite set S, a distribution \i over S is a function fx : S — > [0, 1] satisfying the 
condition fx(S) < 1. The support of fi is defined by Supp(fi) = {s | /i(s) > 0}, and the 
size of n is defined by \n\ = \Supp(fi)\. The distribution p, is called stochastic if fi(S) = 1, 
absorbing if p(S) = 0. We sometimes use an auxiliary state (not a real state) _L g" S and 
set n(-L) = 1 — n(S). If [i is not stochastic we have //(-L) > 0. Further, let S± denote the 
set S U {-L}, and let Supp ± (fi) = Supp(fi) U {_L} if /i(^) > and Supp ± (fi) = Supp(fi) 
otherwise. We let Dist(S) denote the set of distributions over the set S. 

Definition 2.1. A labelled fully probabilistic system (FPS) is a tuple A4 = (S, P, L) where 
S is a finite set of states, P : S x S — > [0, 1] is a probability matrix such that P(s, •) £ Dist(S) 
for all s € S, and L : S — > 2 AP is a labelling function. 

A state s is called stochastic and absorbing if the distribution P(s, •) is stochastic 
and absorbing respectively. For s £ S, let post(s) = Supp(P(s, •)), and let post^(s) = 
Supp ± (P(s, •))• 

Definition 2.2. A labelled discrete-time Markov chain (DTMC) is an FPS M = (S,P,L) 
where s is either absorbing or stochastic for all s £ S. 

FPSs and DTMCs are time- abstract, since the duration between triggering transitions 
is disregarded. We observe the state only at a discrete set of time points 0, 1,2, . . .. We 
recall the definition of CTMCs which are time-aware: 

Definition 2.3. A labelled continuous-time Markov chain (CTMC) is a tuple A4 = (S, R, L) 
with S and L as before, and R:SxS^ M>o is a rate matrix. 

For CTMC M, let post(s) = {s' £ S | R(s,s') > 0} for all s £ S. The rates give 
the average delay of the corresponding transitions. Starting from state s, the probability 
that within time t a successor state is chosen is given by 1 — e - R -( s > s ')*. The probability 
that a specific successor state s' is chosen within time t is thus given by (1 — e -R "( s ' 5 ')') • 
R(s, s')/R(s, S). A CTMC induces an embedded DTMC, which captures the time-abstract 
behaviour of it: 
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Definition 2.4. Let M = (S, R, L) be a CTMC. The embedded DTMC of M is defined 
by emb(M) = (S,P,L) with P(s,s') = R(s, s')/R(s, S) if R(s,S) > and otherwise. 

We will also use P for a CTMC directly, without referring to its embedded DTMC 
explicitly. If one is interested in time-abstract properties (e.g., the probability to reach a 
set of states) of a CTMC, it is sufficient to analyse its embedded DTMC. 

For a given FPS, DTMC or CTMC, its fanout is defined by max s6 5 \post(s)\. The 
number of states is defined by n = \S\, and the number of transitions is defined by m = 
SsgS \post(s)\. For s G S, reach(s) denotes the set of states that are reachable from s 
with positive probability. For a relation R C S x S and s G S, let R(s) denote the set 
{s 1 G S | (s,s') G R}. Similarly, for s' G S, let i?" 1 ^') denote the set {s G 5 | (s,s') G i?}. 
If (s, s') G i?, we write also s R s'. 

Markov chains are purely probabilistic. Now we consider extensions of Markov chains 
with nondeterminism. We first recall the definition of probabilistic automata, which can be 
considered as the simple probabilistic automata with transitions allowing deadlocks in |39j . 

Definition 2.5. A probabilistic automaton (PA) is a tuple M = (S,Act,P,L) where S 
and L are defined as before, Act is a finite set of actions, PCSx Act x Dist(S) is a finite 
set, called the probabilistic transition matrix. 

For (s, a, /x) G P, we use s —> [i as a shorthand notation, and call /i an a-successor distri- 
bution of s. Let Act(s) = {a \ 3fx : s fi} denote the set of actions enabled at s. For s G S 
and a G Act(s), let Steps a (s) = {fi G Dist(S) \ s //} and Steps(s) = UoGAct(s) Steps a (s). 
The fanout of a state s is defined by /<m(s) = £ aG Arf(s) E Me step Sa ( s )(H + 1)- Intuitively, 
fan(s) denotes the total sum of the sizes of outgoing distributions of state s plus their 
labelling. The fanout of A4 is defined by max s gg fan(s). Summing up over all states, we 
define the size of the transitions by m = X^sgs f an ( s )- 

A Markov decision process (MDP) [36] arises from a PA A4 if for s £ S and a G Act, 
there is at most one a-successor distribution \x of s, which must be stochastic. 

We consider a continuous-time counterpart of PAs where the transitions are described 
by rates instead of probabilities. A rate function is simply a function r : S — > M>o- Let 
\r\ = \{s | r(s) > 0}| denote the size of r. Let Rate(S) denote the set of all rate functions. 

Definition 2.6. A continuous-time PA (CPA) is a tuple (S, Act, R, L) where S, Act, L as 
defined for PAs, and RCSx Act x Rate(S) a finite set, called the rate matrix. 

We write s —* r if (s, a, r) G R, and call r an a-successor rate function of s. For 
transition s — > r, the sum r(S) is also called the exit rate of it. Given that the transition 
s —> r is chosen from state s, the probability that any successor state is chosen within time 
t is given by 1 — e~ r ^ t , and a specific successor state s' is chosen within time t is given by 
(1 — e~ r ^ 1 ) ■ T^y- The notion of Act(s), Steps a {s), Steps(s), fanout and size of transitions 
for PAs can be extended to CPAs by replacing occurrence of distribution \i by rate function 
r in an obvious way. 

The model continuous-time Markov decision processes (CTMDPs) [361 [7] can be con- 
sidered as special CPAs where for s G S and a G Act, there exists at most one rate function 
r G Rate(S) such that s — > r. The model CTMDPs considered in paper [32] essentially 
agree with our CPAs. 
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Figure 1: An FPS for illustrating the simulation relation^. 



2.2. Strong and Weak Simulation Relations. We first recall the notion of strong sim- 
ulation on Markov chains [9], PAs [ID], and CPAs [Mj. Strong probabilistic simulation is 
defined in Subsection 12.2.21 Weak simulation for Markov chains will be given in Subsec- 
tion 12.2.31 The notion of simulation up to R is introduced in Subsection 12.2.41 

2.2.1. Strong Simulation. Strong simulation requires that each successor distribution of one 
state has a corresponding successor distribution of the other state. The correspondence of 
distributions is naturally defined with the concept of weight functions [26j , adapted to FPSs 
as in [9j. For a relation R C S x S, we let R± denote the set R U {(-L, s) \ s G S±}. 

Definition 2.7. Let fj,, fjf € Dist(S) and R C S x S. A weight function for (/i,/U ; ) with 
respect to R is a function A : S± x S± — > [0, 1] such that 

(1) A(s, s') > implies s R± s', 

(2) n(s) = A(s,S±) for s € S± and 

(3) /i'(s') = A(S ± ,s') for s' £ 5±. 

We write [i Qr fi' if there exists a weight function for (//, //) with respect to R. 

The first condition requires that only pairs (s,s') in the relation R± have a positive 
weight. In other words, for s, s' € 5 with s' R±(s), it holds that A(s, s') = 0. Strong sim- 
ulation requires similar states to be related via weight functions on their distributions [26] . 

Definition 2.8. Let M = (S,P,L) be an FPS, and let R C S x S. The relation i? is a 
strong simulation on A4 iff for all si, S2 with si S2: ^(si) = L{s2) and P(sj, •) Qr P(s2, •)• 
We say that S2 strongly simulates si in A4. denoted by s± s 2, iff there exists a 
strong simulation R on M such that si S2- 

By definition, it can be shown [9J that is reflexive and transitive, thus a preorder. 
Moreover, is the coarsest strong simulation relation for Ai. If the model A4 is clear from 
the context, the subscript M. may be omitted. Assume that s\ ^ S2 and let A denote the 
corresponding weight function. If P(s2, -L) > 0, we have that P(s2, -L) = J2 s eS ^( s ' -L) = 
A(_L, _L). The second equality follows by the fact that _L can not strongly simulate any real 
state in S. Another observation is that if s is absorbing, then it can be strongly simulated 
by any other state s' with L(s) = L(s'). 



Although this graph is not connected, it shows a single FPS. Similarly, later figures will show a single 
DTMC, CTMC etc. 
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Example 2.9. Consider the FPS depicted in Figure [H Recall that labelling of states is 
indicated by colours in the states. Since the yellow (grey) states are absorbing, they strongly 
simulate each other. The same holds for the green (dark grey) states. We show now that 
si 3 s 2 but s 2 £ s 3 . 

Consider first the pair (si,S2). Let R = {(si, S2), (1*1, u 2 ), {v±, V2), (qi, 92)}- We show 
that R is a strong simulation relation. First observe that L(s) = L(s') for all (s, s') € R. 
Since states u\,q\ are absorbing, the conditions for the pairs (141,1*2) and (q\,q 2 ) hold 
trivially. To show the conditions for (v\, V2), we consider the function Ai defined by: 
Ax(gi,q , 2) = gj Ai(J_,g3) = |, Ai(_L,_L) = g and Ai(-) = otherwise. It is easy to check 
that Ai is a weight function for (P(ui, •), P(i>2, ■)) with respect to R. Now consider (si, S2). 
The weight function A 2 for (P(si, •),P(s2 ) ')) with respect to it! is given by A2 (1*1,1x2) = \ 
and /S.2{vi,V2) = A.2(-L,V2) = \ and A 2 {-) = otherwise. Thus R is a strong simulation 
which implies that s\ ^ s 2 - 

Consider the pair (s 2 ,ss). Since P(s2,«2) = 2' to establish the Condition [2] of Defini- 
tion [221 we should have \ = A(v 2 ,S±). Observe that V3 is the only successor state of S3 
which can strongly simulate i>2, thus A(i)2,S±) = A(w2>i>3)- However, for state V3 we have 
P(s3,i'3) < A(v2,vs), which violates the Condition l3l of Definition 12.71 thus we cannot find 
such a weight function. Hence, S2 j£ S3. 

Since each DTMC is a special case of an FPS, Definition 12 .81 applies directly for DTMCs. 
For CTMCs we say that S2 strongly simulates s\ if, in addition to the DTMC conditions, 
S2 can move stochastically faster than s\ [9], which manifests itself by a higher rate. 

Definition 2.10. Let M = (S,K,L) be a CTMC and let R C S x S. The relation R is a 
strong simulation on M. iff for all si,S2 with si R S2'- L(s\) = L(s2), P(si,-) P(s2,-) 
and R(si,5) < H(s 2 ,S). 

We say that S2 strongly simulates si in Ai, denoted by si S2, iff there exists a 
strong simulation R on M. such that si R s 2 - 

Thus, s s' holds if s ^ e mb(M) s 'i an d s> 1S faster than s. By definition, it can be 
shown that is a preorder, and is the coarsest strong simulation relation for A4. For 
PAs, strong simulation requires that every a-successor distribution of s\ is related to an 
a-successor distribution of S2 via a weight function \4Q\ [26] : 

Definition 2.11. Let M = (S, Act, P, L) be a PA and let R C S x S. The relation R is a 
strong simulation on M iff for all si,S2 with si R S2- L(s\) = L{s2) and if si —* [i\ then 
there exists a transition s 2 /-*2 with /ii C^j /i 2 . 

We say that S2 strongly simulates s% in AI, denoted si s 2, iff there exists a strong 
simulation R on A4 such that si R S2- 

Example 2.12. Consider the PA in Figure [2[ Then, it is easy to check si ^ S2- each 
a-successor distribution of s± has a corresponding a-successor distribution of S2- However, 
si does not strongly simulate s 2 , as the middle a-successor distribution of S2 can not be 
related by any a-successor distribution of s±. 

Now we consider CPAs. For a rate function r, we let /i(r) € Dist(S) denote the 
induced distribution defined by: if r(S) > then /j(r)(s) equals r(s)/r(S) for all s £ S, and 
if r(S) = 0, then **(r)(s) = for all s € S. Now we introduce the notion of strong simulation 
for CPAs [H], which can be considered as an extension of the definition for CTMCs [9]: 
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Figure 2: A PA for illustrating the simulation relations. 



Definition 2.13. Let M = (S, Act, R, L) be a CPA and let R C S x S. The relation R is 
a strong simulation on M. iff for all s±, S2 with s\ R S2- L(si) = L(s2) and if s% —> r% then 
there exists a transition S2 — > ^2 with //(rr) ^(^2) and ?r(5) < ^(S). 

We write si ^» S2 iff there exists a strong simulation R on A4 such that si i? S2. 

Similar to CTMCs, the additional rate condition ri(S) < ^(S 1 ) indicates that the 
transition S2 — ► ^2 is faster than s\ —> 7*1. As a shorthand notation, we use ri r2 for the 
condition /x(r*i) C.r ^(^2) and r\(S) < ^(S*). For both PAs and CPAs, is the coarsest 
strong simulation relation. 

2.2.2. Strong Probabilistic Simulations. We recall the definition of strong probabilistic sim- 
ulation, which is coarser than strong simulation, but still preserves the same class of PCTL- 
properties as strong simulation does. We first recall the notion of combined transition [39], 
a convex combination of several equally labelled transitions: 

Definition 2.14. Let M = (S,Act,P,L) be a PA. Let s G S, a G Act(s) and k = 
\Steps Q (s)\. Assume that Steps a (s) = {//!,...,//&}. The tuple (s,a,/i) is a combined 
transition, denoted by s ~* fi, iff there exist constants ci, . . . , G [0, 1] with X^i=i c « = 1 
such that /x = Yli=i CifH- 

The key difference to Definition 12.111 is the use of ~* instead of 

Definition 2.15. Let M = (S,Act,P,L) be a PA and let R C S 1 x 5. The relation i? is 
a strong probabilistic simulation on A4 iff for all si,S2 with si R S2- L(si) = L{s2) and if 
si — > m then there exists a combined transition S2 ^ ^2 with //1 C fi 

We write s\ S2 iff there exists a strong probabilistic simulation R on A4 such that 
si it! s 2 - 

Strong probabilistic simulation is insensitive to combined transitions^, thus, it is a 
relaxation of strong simulation. Similar to strong simulation, is the coarsest strong 
probabilistic simulation relation for A4. Since MDPs can be considered as special PAs, 
we obtain the notions of strong simulation and strong probabilistic simulation for MDPs. 
Moreover, these two relations coincide for MDPs as, by definition, for each state there is at 
most one successor distribution per action. 

The combined transition denned in [39] is more general in two dimensions: First, successor distributions 
are allowed to combine different actions. Second, Yl =i c» < 1 is possible. The induced strong probabilistic 
probabilistic simulation preorder is, however, the same. 
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Figure 3: A Continuous-time Probabilistic Automaton. 



Example 2.16. We consider again the PA depicted in Figure [21 From Example 12.121 we 
know that S2 ^ s\. In comparison to state si, state S2 has one additional a-successor 
distribution: to states and with equal probability 0.5. This successor distribution can 
be considered as a combined transition of the two successor distributions of si: each with 
constant 0.5. Hence, we have S2 ^ p s\. 

We extend the notion of strong probabilistic simulation for PAs to CPAs. First, we 
introduce the notion of combined transitions for CPAs. In CPAs the probability that a 
transition occurs is exponentially distributed. The combined transition should also be 
exponentially distributed. The following example shows that a straightforward extension of 
Definition 12.141 does not work. 

Example 2.17. For this purpose we consider the CPA in Figure [3l Let r\ and ri denote 
left and the right a-successor rate functions out of state s\. Obviously, they have different 
exit rates: r\{S) = 10, r2(S) = 18. Taking each with probability 0.5, we would get the 
combined transition r = 0.5ri + 0.5r2: r({ui,U2}) = 7 and r({vi,V2}) = 7. However, 
r is hyper-exponentially distributed: the probability of reaching yellow (grey) states (ui 
or U2) within time t under r is given by: 0.5 • ^ • (1 — e~ 10< ) + 0.5 • y| • (1 — e -18 *). 
Similarly, the probability of reaching green (dark grey) states within time t is given by: 
0.5 • ^ • (1 - e" 10 *) + 0.5 • ^ • (1 - e" 18 *). 

From state S2, the two a-successor rate functions have the same exit rate 14. Let r[ 
and r' 2 denote left and the right a-successor rate functions out of state S2- In this case the 
combined transition r' = 0.5r[ + 0.5r 2 is also exponentially distributed with rate 14: the 
probability to reach yellow (grey) states («3 and 114) within time t is X • (1 — e -14 '), which 
is the same as the probability of reaching green (dark grey) states within time t. 

Based on the above example, it is easy to see that to get a combined transition which 
is still exponentially distributed, we must consider rate functions with the same exit rate: 

Definition 2.18. Let M = (S,Act,K,L) be a CPA. Let s £ S, a € Act(s) and let 
{ri,...,rfc} C Steps a (s) where ri(S) = rj(S) for i,j £ {l,...,k}. The tuple (s,a,r) is 
a combined transition, denoted by s "2* r, iff there exist constants c\, . . . ,Cf- G [0, 1] with 

Yji=1 °i = 1 sucn tllat r = S»=l °i r i- 

In the above definition, unlike for the PA case, only a-successor rate functions with 
the same exit rate are combined together. Similar to PAs, strong probabilistic simulation 
is insensitive to combined transitions, which is thus a relaxation of strong simulation: 

Definition 2.19. Let M = (S, Act, R, L) be a CPA and let R C S x S. The relation R is 
a strong probabilistic simulation on A4 iff for all si,S2 with s± R S2- L(si) = L{s2) and if 
s\ T\ then there exists a combined transition S2 ~> T2 with T\ Qr T2- 
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Figure 4: Splitting of successor states for weak simulations. 



We write si S2 iff there exists a strong simulation R on Ai such that s\ R S2- 

Recall 7*1 Qr T2 is a shorthand notation for \x{v\) Qr /^(r^) and r\{S) < r2(S). By 
definition, the defined strong probabilistic simulation ~^ P M is the coarsest strong probabilistic 
simulation relation for Ai. 

Example 2.20. Reconsider the CPA in Figure [3l As discussed in Example 12.171 the two 

a-successor rate functions of s\ cannot be combined together, thus the relation sq ^ p s\ 
cannot be established. However, sq ^ p S2 holds: denoting the left rate function of S2 as r\ 
and the right rate function as r2, we choose as the combined rate function r = 0.5ri + 0.5r2. 
Obviously, the conditions in Definition 12.191 are satisfied. 

2.2.3. Weak Simulations. We now recall the notion of weak simulation [9] on Markov 
chains^. Intuitively, S2 weakly simulates s\ if they have the same labelling, and if their 
successor states can be grouped into sets Ui and V{ for i = 1,2, satisfying certain condi- 
tions. Consider Figure HI We can view steps to Vj as stutter steps while steps to Ui are 
visible steps. With respect to the visible steps, it is then required that there exists a weight 
function for the conditional distributions: P(si, -)/Ki and P(s2, ")/-^2 where ifj intuitively 
is the probability to perform a visible step from Sj. The stutter steps must respect the 
weak simulation relations: thus states in V2 should weakly simulate s\, and state S2 should 
weakly simulate states in V\. This is depicted by dashed arrows in the figure. For reasons 
we will explain later in Example 12.22^ the definition needs to account for states which par- 
tially belong to Ui and partially to Vi- Technically, this is achieved by functions Si that 
distribute Sj over Ui and V% in the definition below. For a given pair (s\,S2) an d functions 
Si : S — > [0, 1], let Us^ V$ z C S (for i = 1,2) denote the sets 

U Si = {u£ post{ Si ) I Si{u) > 0}, V Si = {v £ post(si) I Si(v) < 1} (2.1) 

If (si, S2) and S( are clear from the context, we write Ui, Vi instead. 

Definition 2.21. Let M = (S,P,L) be a DTMC and let R C S x S. The relation R is a 
weak simulation on M. iff for all s\, S2 with s\ R S2- L(s\) = L{s2) and there exist functions 
Si : S — > [0, 1] such that: 

(1) (a) v\ R S2 for all v\ E V\, and (b) S\ R V2 for all V2 £ V2 

(2) there exists a function A : S x S — > [0, 1] such that: 

(a) A(ui,U2) > implies u\ G U\,U2 S U2 and ui R 1*2- 

4 In [45], we have also considered decision algorithm for weak simulation for FPSs, which is defined in [9]. 
However, as indicated in [33], the proposed weak simulation for FPSs contains a subtle flaw, which cannot 
be fixed in an obvious way. Thus, in this paper we restrict to weak simulation on DTMCs and CTMCs. 
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Figure 5: A DTMC where splitting states is necessary to establish the weak simulation. In 
the model some states are drawn more than once. 

(b) if K\ > and K 2 > then for all states w G S: 

K\ ■ A(w,U 2 ) = P(s 1 ,w)5 1 (w), K 2 ■ A(Ui,w) = P(s 2 ,w)5 2 (w) 

where K,, = Y J u i au x ?>i( u i) ' P{ s i, u i) f° r * = 1, 2. 
(3) for u\ € U\ there exists a path fragment s 2 , W\, . . . , w ni u 2 with positive probability such 

that n > 0, si R Wj for < j < n, and u\ R u 2 . 
We say that s 2 weakly simulates s% in A4, denoted si T^m s 2 , iff there exists a weak 
simulation R on M. such that si i? s 2 . 

Note again that the sets Ui,Vi in the above definition are defined according to Equa- 
tion [2J] with respect to the pair (si,s 2 ) and the functions <5j. The functions 6{ can be 
considered as a generalisation of the characteristic function of Ui in the sense that we may 
split the membership of a state to Ui and Vi into fragments which sum up to 1. For exam- 
ple, if Si(s) = |, we say that 3 fragment of the state s belongs to Ui, and | fragment of s 
belongs to V\. Hence, Ui and Vi are not necessarily disjoint. Observe that Ui = implies 
that 5i(s) = for all s € 5. Similarly, = implies that <5j(s) = 1 for all s £ S. 

Condition [3] will in the sequel be called the reachability condition. If K% > and 
K 2 = 0, which implies that U 2 = and ?7i ^ 0, the reachability condition guarantees that 
for any visible step si — > u% with u\ G Ui, s 2 can reach a state u 2 which simulates ui while 
passing only through states simulating s\. Assume that we have S = {si, s 2 ,ui} where 
L(si) = L{s 2 ) and u\ has a different labelling. There is only one transition P(si,ui) = 1. 
Obviously s\ ^ s 2 . Dropping Condition [3] would mean that si ^ s 2 . We illustrate the use 
of fragments of states in the following example: 

Example 2.22. Consider the DTMC depicted in Figure [5j For states m, u 2 , vi, v 2 , obvi- 
ously the following pairs (ui,u 2 ), (ui,v 2 ), (vi,v 2 ) are in the weak simulation relation. The 
state u 2 cannot weakly simulate vi. Since v 2 weakly simulates vi, it holds that s 2 ^ S5. 
Similarly, from ui ^ vi we can easily show that si ^ S4. We observe also that s 2 ^ S3: 
Ki > and K 2 > since both s 2 and S3 have yellow (grey) successor states, but the 
required function A cannot be established since u 2 cannot weakly simulate any successor 
state of s 2 (which is vi). Thus s 2 ^ S3. 
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Without considering fragments of states, we show that a weak simulation between si 
and S3 cannot be established. Since S2 ^ S3, we must have U\ = {u\,vi,S2} and V\ = 0. 
The function 5\ is thus defined by 5\{u\) = 8\{v\) = <5(s2) = 1 which implies that K\ = 1. 
Now consider the successor states of S3. Obviously #2(^2) = ^2(^2) = 1> which implies that 
U2,V2 G U2. We consider the following two cases: 

• The case 52(54) = 1- In this case we have that K2 = 1. A function A must be defined 
satisfying Condition I2bl in Definition 12.211 Taking w = S4, the following must hold: 
K 2 ■ A(Ui,s 4 ) = P(s 3 ,s 4 )J 2 (s 4 ). As K 2 = l,P(s 3 ,s 4 ) = 0.75 and 5 2 (s 4 ) = 1, it follows 
that A(C/i,s 4 ) = 0.75. The state S2 is the only successor of s 4 that can be weakly 
simulated by s 4 , so A(s2,s 4 ) = 0.75 must hold. However, the equation K\ ■ A(s2,U2) = 
P(si) •S2)^i('S2) does not hold any more, as on the left side we have 0.75 but on the right 
side we have 0.5 instead. 

• The case #2(s 4 ) = 0. In this case we have still K2 > 0. Similar to the previous case it is 
easy to see that the required function A cannot be defined: the equation K\ ■ A(s2, U2) = 
P( s i) S2)5\(s2) does not hold since the left side is (no states in U2 can weakly simulate 
S2) but the right side equals 0.5. 

Thus without splitting, S3 does not weakly simulate s 4 . We show it holds that s 4 ^ S3. It is 
sufficient to show that the relation R = {(si, S3), (u\ } U2), (i>i, 1*2)5 (q±, ^3), (s 4 , s 4 ), (u\,v\), 
(vx,vi),(qi,qi), (S2, S5), (S2, s 4 )} is a weak simulation relation. By the discussions above, 
it is easy to verify that every pair except (si, S3) satisfies the conditions in Definition 12.211 
We show now that the conditions hold also for the pair (s\,ss). The function 5\ with 
5\{u\) = 5i{v\) = 5i(s2) = 1 is defined as above, also the sets U\ = {1*1,^1,52}, ^1 = 
The function 82 is defined by: #2(^2) = #2(^2) = 1 and #2(s 4 ) = |, which implies that 
U2 = {u2,V2,Si} and V2 = {s 4 }. Thus, we have K\ = 1 and K2 = 0.5. Since s\ ^ s 4 , 
Condition [1] holds trivially as (s 4 ,s 4 ) G R. The reachability condition also holds trivially. 
To show that Condition [2] holds, we define the function A as follows: A(ui,U2) = 0.3, 
A(ui,u 2 ) = 0.2 and A(s 2 ,s 4 ) = 0.5. We show that K 2 ■ A(Pi,u;) = P(s 3 , w)S 2 (w) holds 
for all w G S. It holds that K2 = 0.5. First observe that for w G" U2 both sides of the 
equation equal 0. Let first w = U2 for which we have that P(s3, ^2)52(^2) = 0.15. Since 
A(U\,U2) = 0.3, also the left side equals 0.15. The case w = V2 can be shown in a similar 
way. Now consider w = s 4 . Observe that A(£/i, s 4 ) = 0.5 thus the left side equals 0.25. The 
right side equals P(s3, s 4 )<52(s 4 ) = 0.75 • | = 0.25 thus the equation holds. The equation 
K\ ■ A(w,U2) = P(si, w)5i(w) can be shown in a similar way. Thus A satisfies all the 
conditions, which implies that si ^ S3. 

Weak simulation for CTMCs is defined as follows. 

Definition 2.23 (PE]). Let M = (S,K,L) be a CTMC and let R C S x S. The relation 
R is a weak simulation on A4 iff for s 4 R 82- L(si) = L(s2) and there exist functions 
5i : S — > [0, 1] (for i = 1, 2) satisfying Equation 12. II and Conditions [T] and [2] of Definition 12.2 II 
and the rate condition: 

(3') Kx-R( Sl ,S) <K 2 -R(s 2 ,S) 

We say that S2 weakly simulates s 4 in A4, denoted si s 2, iff there exists a weak 

simulation R on A4 such that s 4 R S2- 

In this definition, the rate condition [3j strengthens the reachability condition of the 
preceding definition. If U\ ^ 0, we have that K\ > 0; the rate condition then requires that 
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Figure 6: A simple FPS for illustrating the simulation up to R. 



K 2 > 0, which implies U 2 / 0. For both DTMCs and CTMCs, the defined weak simulation 
^ is a preorder [9], and is the coarsest weak simulation relation for A4. 

2.2.4. Simulation up to R. For an arbitrary relation R on the state space S of an FPS with 
s\ R s 2 , we say that s 2 simulates s± strongly up to R, denoted s\ ^r s 2 , if L(si) = L(s 2 ) 
and P(si, •) Qr P(s 2 , •)■ Otherwise we write s± s 2 . Since only the first step is considered 
f° r ^fii s i ^-R s 2 does not imply s± s 2 unless R is a strong simulation. By definition, 
R is a strong simulation if and only if for all s\ R s 2 it holds that s\ ^r s 2 . Likewise, we 
say that s 2 simulates s\ weakly up to R, denoted by si ^r s 2 , if there are functions 5{ and 
Ui,Vi,A as required by Definition 1 2 . 2 1 1 for this pair of states. Otherwise, we write s\ ^£r s 2 . 
Similar to strong simulation up to R, s\ ^.r s 2 does not imply s\ s 2 , since no conditions 
are imposed on pairs in R different from (si, £2). Again, R is a weak simulation if and only 
if for all s\ R s 2 it holds that s\ t6r s 2 . These conventions extend to DTMCs, CTMCs, 
PAs and CPAs in an obvious way. For PAs and CPAs, strong probabilistic simulation up 
to R, denoted by ~^ R , is also defined analogously. 

Example 2.24. Consider the FPS in Figured Let R = {(si,s 2 ), (w\,w 2 )}. Since L(qi) ^ 
L(q 2 ) we have that w\ w 2 . Thus, R is not a strong simulation. However, s± ^r s 2 , as 
the weight function is given by A(wi,w 2 ) = 1. Let R' = {(si,s 2 )}, then, s± ^ri s 2 . 



3. Maximum Flow Problems 

Before introducing algorithms to decide the simulation preorder, we briefly recall the preflow 
algorithm [20j for finding the maximum flow over the network M = (V, E, c) where V is a 
finite set of vertices, E C V x V is a set of edges, and c : E — > M>o U {00} is the capacity 
function. V contains a distinguished source vertex / and a distinguished sink vertex \. We 
extend the capacity function to all vertex pairs: c(v, w) = if (v, w) E. A flow f on M 
is a function / : V x V — > M that satisfies: 

(1) f(v,w) < c(v,w) for all (v, w) G V x V capacity constraints 

(2) f(v,w) = —f(w,v) for all (v,w) G V x V antisymmetry constraint 

(3) f(V, v) = at vertices v £ V \ {f~ ,\} conservation rule 
The value of a flow function / is given by /(/, V). A maximum flow is a flow of maximum 
value. A preflow is a function / : V x V — > M satisfying Conditions [Q and [2] above, and the 
relaxation of Condition [3j 

(13) /(V, v) > for all v G V \ {/}. 

The excess e(v) of a vertex v is defined by f(V,v). A vertex v {/A} is called active if 
e(u ) > 0. Observe that if in a preflow function no vertex v is active for v G V \ {/*, \}, it 
is then also a flow function. A pair (v,w) is a residual edge of / if f(v,w) < c(v,w). The 
set of residual edges with respect to / is denoted by Ef. The residual capacity cf(v,w) of 
the residual edge (v,w) is defined by c(v,w) — f(v,w). If (v,w) is not a residual edge, it 
is called saturated. A valid distance function (called valid labelling in [20]) d is a function 
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V — > N U {00} satisfying: d(f) = \V\, d(\) = and d(v) < d(w) + 1 for every residual edge 
(v, w). A residual edge (v, w) is admissible if d(v) = d(w) + 1. 

Related to maximum flows are minimum cuts. A cut of a network J\f = (V, E, c) is a 
partition of V into two disjoint sets (X, X') such that / € X and \€l'. The capacity of 
(X, X') is the sum of all capacities of edges from X to X', i. e., J2 v ex wex' c ( u ' ^ ™ n " 
imum cut is a cut with minimal capacity. The Maximum Flow Minimum Cut Theorem pQ 
states that the capacity of a minimum cut is equal to the value of a maximum flow. 

The Preflow Algorithm. We initialise the preflow / by: f(v,w) = c(v,w) if v = f and 
otherwise. The distance function d is initialised by: d(v) = \V\ if v = / and otherwise. 
The preflow algorithm preserves the validity of the preflow / and the distance function 
d. If there is an active vertex v such that the residual edge (v, w) is admissible, we push 
5 := min{e(f ), Cf{v , w)} amount of flow from v toward the sink along the admissible edge 
(y, w) by increasing f(v, w) (and decreasing f(w, v)) by 5. The excesses of v and w are then 
modified accordingly by: e(v) = e(v) — 5 and e(w) = e{w) + 5. If v is active but there are no 
admissible edges leaving it, one may relabel v by letting d(v) := mm{d(w) + l | (v, w) G Ef}. 
Pushing and relabelling are repeated until there are no active vertices left. The algorithm 
terminates if no such operations apply. The resulting final preflow / is a maximum flow. 

Feasible Flow Problem. Let A C E be a subset of edges of the network M = (V, E, c), and 
define the lower bound function I : A — > M>o which satisfies 1(e) < c(e) for all e E A. We 
address the feasible flow problem which consists of finding a flow function / satisfying the 
condition: /(e) > 1(e) for all e E A. We briefly show that this problem can be reduced to 
the maximum flow problem pQ. 

We can replace a minimum flow requirement on edge v — > w by turning v into a 
demanding vertex (i.e., a vertex that consumes part of its inflow) and turning w into a 
supplying vertex (i.e., a vertex that creates some outflow ex nihilo). The capacity of edge 
v — > w is then reduced accordingly. 

Now, we are going to look for a flow-like function for the updated network. The function 
should satisfy the capacity constraints, and the difference between outflow and inflow in 
each vertex corresponds to its supply or demand, except for / and V To remove that last 
exception, we add an edge from \ to / with capacity 00. 

We then apply another transformation to the updated network so that we can apply 
the maximum flow algorithm. We add new source and target vertices /' and \' . For each 
supplying vertex s, we add an edge /' — > s with the same capacity as the supply of the 
vertex. For each demanding vertex d, we add an edge d — ► V with the same capacity as 
the demand of the vertex. In [lj it is shown that the original network has a feasible flow if 
and only if the transformed network has a flow h that saturates all edges from /' and all 
edges to \, . The flow h necessarily is a maximum flow, and if there is an h, each maximum 
flow satisfies the requirement; therefore it can be found by the maximum flow algorithm. 
An example will be given in Example 15.41 in Section [5l 

4. Algorithms for Deciding Strong Simulations 

We first recall the basic algorithm to compute the largest strong simulation relation ^ in 
Subsection 14.11 Then, we refine this algorithm to deal with strong simulation on Markov 
chains in Subsection l4.21 and extend it to deal with probabilistic automata in Subsection l4.31 
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SimRel s (A1) 

1.1: Rx <- {(si, s 2 ) £ S x S \ L(si) = L(s 2 )} and i <- 

1.2: repeat 

1.3: i i + 1 

1.4: <- 

1.5: for all (si,S2) £ do 

1.6: if si s 2 then 

1.7: fl i+ i<-i2 i+ iU{(ai,a 2 )} 
1.8: until Ri+i = Ri 
1.9: return i?. 

Algorithm 1: Basic algorithm to decide strong simulation. 

In Subsection 14.41 we present an algorithm for deciding strong probabilistic simulation for 
probabilistic automata. 

4.1. Basic Algorithm to Decide Strong Simulation. The algorithm in [3j, copied as 
SimReLs in Algorithm [H takes as a parameter a model, which, for now, is an FPS Ai. The 
subscript 's' stands for strong simulation; a very similar algorithm, namely SimRel^, will 
be used for weak simulation later. To calculate the strong simulation relation for Ai, the 
algorithm starts with the initial relation R\ = {(si,«2) € S X S \ L(si) = L(s2)} which is 
coarser than ^m- m iteration i, it generates Ri+i from Ri by deleting each pair (si,s 2 ) 
from Ri if S2 cannot strongly simulate s\ up to Ri, i.e., si S2- This proceeds until 
there is no such pair left, i.e., Ri+\ = Ri. Invariantly throughout the loop it holds that 
Ri is coarser than (i.e., is a sub-relation of Ri). We obtain the strong simulation 
preorder = Ri, once the algorithm terminates. 

The decisive part of the algorithm is the check in Line ll.6( i. e., whether s% 82- This 
can be answered via solving a maximum flow problem on a particular network A/"(P(si, •), 
P(s2, -),Ri) constructed from P(si, ■), P(s2, •) and Ri. This network is the relevant part of 
a graph containing two copies t G S± and t G S± of each state where S± = {t \ t G S±} 
as follows: Let / (the source) and \ (the sink) be two additional vertices not contained in 
S±L)S±. For [i, fjf £ Dist(S), and a relation R C S x S we define the network M(fi, fx', R) = 
(V,E,c) with the set of vertices V = {/A} U Supp^(/j,) U Suppj_(/i') and the set of edges 
E is defined by E = {(s,t) \ (s,t) eR±} U {(As)} U {(t,\)} where s G Supp ± (fi) and 
t G Suppj_(fi'). Recall the relation R± is defined by R U {(_L,s) | s G The capacity 

function c is defined as follows: c(/,s) = ;u(s) for all s G Suppj_(fi), c(t,\) = n'(t) for 
all t G Supp^fi'), and c(s,f) = 00 for all other (s,i) G £7. This network is a bipartite 
network, since the vertices can be partitioned into two subsets V\ := Supp^(n) U {\} 
and V2 := Suppj_([i') U {/} such that all edges have one endpoint in V\ and another in 
V2. Later, we will use two variations of this network: For 7 G K>o, we let N{n, 7//', R) 
denote the network obtained from Af(fi, fi' , R) by setting the capacities to the sink \ to: 
c(t, \) = 7/x'(t). For two states s\, S2 of an FPS or a CTMC, we let M(s±, S2, R) denote the 
network AA(P(si, •), P(s2, -),R)- 

The following lemma expresses the crucial relationship between maximum flows and 
weight functions on which the algorithm is based. It is a direct extension of [3j Lemma 5.1]: 
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Lemma 4.1. Let S be a finite set of states and R be a relation on S. Let G Dist(S). 
Then, \x fj, 1 iff the maximum flow of the network Af(/J,, fi', R) has value 1. 

Proof. As we introduced the auxiliary state _L, /U and are stochastic distributions in 
Dist(S±). The rest of the proof follows directly from [3j Lemma 5.1]. □ 

Thus we can decide s\ 7^ S2 by computing the maximum flow in M(si, S2, Ri) and 
then check whether it has value 1. We recall the correctness and complexity of SimRel s 
which will also be used later. 

Theorem 4.2 ([3]). // SimRel s (A4) terminates, the returned relation equals ^m- More- 
over, SimRel s (A / () runs in time 0{n /logn) and in space 0(n 2 ). 

Proof. First we show that after the last iteration (say iteration A;), it holds that 77; is coarser 
than Rk'. It holds that Rk+i = Rk, thus for all (51,^2) £ Rk, we have that s\ S2- As 
for all (si,S2) £ Rk f= -^l) we have L(s\) = L(s2), Rk is a strong simulation relation by 
Definition 12.81 thus 70 is coarser than Rk- 

Now we show by induction that the loop of the algorithm invariantly ensures that Ri 
is coarser than Assume i = 1. By definition of strong simulation, s± 70 S2 implies 
L(si) = L{s2). Thus, the initial relation R\ is coarser than the simulation relation 7^. Now 
assume that Ri is coarser than 77; for some 1 < i < k; we will show that also Ri+i is coarser 
than 73. Pick a pair (si, S2) € 70 arbitrarily. By Definition 12.81 P(si, •) C-< P(s2> ■), so there 
exists a weight function for (P(si, •), P(s2, •)) with respect to 7> Inspection of Definition 12. 71 
shows that the same function is also a weight function with respect to any set coarser than 
7> As Ri is coarser than 70 by induction hypothesis, we conclude that P(si, ■) Qn t P(s2, •), 
and from Subsection 12.2.41 s± 7^ S2- This implies that (si,S2) G Ri+i by line 11.61 for all 
s\ 7^ S2- Therefore, Ri+i is coarser than 70 for all* = 1 . . . , k. 

Now we show the complexity. For one network M(si, S2, Ri) = (V,E,c), the sizes of 
the vertices \V\ and edges \E\ are bounded by 2n + 4 and (n + l) 2 + 2n, respectively. The 
number of edges meets the worst case bound 0(n 2 ). To the best of our knowledge, the best 
complexity of the flow computation for the network G is 0{\ V\ 3 / log \ V\) = 0(n 3 / log n) |14t 
119] . In the algorithm SimRel s , only one pair, in the worst case, is removed from Ri in 
iteration i, which indicates that the test whether si 7^ S2 is called |i?i| times, — 1 

times and so on. Altogether it is bounded by X^i=i * — Y17=i * ^ Oin 4 ). Hence, the overall 
time complexity amounts to C(n 7 /logn). The space complexity is 0(n 2 ) because of the 
representation of the transitions in N(s\, S2, Ri). □ 

4.2. An Improved Algorithm for FPSs. We first analyse the behaviour of SimRel s 
in more detail. For this, we consider an arbitrary pair (s\,S2), and assume that (si,S2) 
stays in relation . . . , Rk throughout the iterations i = 1, . . . , k, until the pair is either 
found not to satisfy s\ S2 or the algorithm terminates with a fix-point after iteration 
k. Then altogether the maximum flow algorithms are run fc-times for this pair. However, 
the networks N(si, S2, Ri) constructed in successive iterations are very similar, and may 
often be identical across iterations: They differ from iteration to iteration only by deletion 
of some edges induced by the successive cleanup of Ri. For our particular pair (51,52) the 
network might not change at all in some iterations, because the deletions from Ri do not 
affect their direct successors. We are going to exploit this observation by an algorithm that 
reuses the already computed maximum flow, in a way that whatever happens is good: If no 
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SMF(i,Af{s 1 ,s 2 ,R i -i), fi-i,di-i, A-i) 

2.1: J\f(si,s 2 ,Ri) <— N(si,S2,Ri-i \ A-l) and /i <- an d <- 

2.2: for all (ui,n 2 ) G A-l do 

2-3: /i(«2,\)^/ife\)-/i(«l,«2) 

2.4: /i(«l,Ba) <- 

2.5: Apply the preflow algorithm to calculate the maximum flow for J\f(s\, s 2 , Ri), 

but initialise the preflow to fi and the distance function to di. 
2.6: return (|/j| = 1, N{s\, s 2 , Ri), fi, di) 

SMF in u{i,S 1 ,S 2 ,R i ) 

2.11: Initialise the network M(si, s 2 , Ri). 

2.12: Apply the preflow algorithm to calculate the maximum flow for M(si, s 2 ,Ri). 
2.13: return = l,Af(s 1 ,s 2 ,R i ), fc,di) 

Algorithm 2: Algorithm for a sequence of maximum flows. 

changes occur from M(si, s 2 , Ri-i) to M(s\, s 2 , Ri), then the maximum flow is the same as 
the one in the previous iteration. If changes do occur, the preflow algorithm can be applied 
to get the new maximum flow very fast, using the maximum flow and distance function 
constructed in the previous iteration as a starting point. 

To understand the algorithm, we look at the network M(si, s 2 , Ri). Let D±, . . . ,Dk 
be pairwise disjoint subsets of R\, which correspond to the pairs deleted from R\ in it- 
eration i, so Ri+i = Ri \ Di for 1 < i < k. Let fj- sl ' S2 ^ denote the maximum flow of 
the network Af(si, s 2 , Ri) for 1 < i < k. We sometimes omit the superscript (si,s 2 ) 
in the parameters if the pair (s\,s 2 ) is clear from the context. We address the prob- 
lem of checking = 1 for all i = 1, . . . , k. Our algorithm sequence of maximum flows 
Smf(z, A/"(si, s 2 , Ri-\), fi-i, di-i, A-l) is shown as Algorithm [2j It executes iteration i of 
a parametric flow algorithm, where A/"(si, s 2 , Ri-i) is the network for (s\, s 2 ) and and 
di-i are the flow and the distance function resulting from the previous iteration % — 1; and 
A-l is a set of edges that have to be deleted from Af(s\, s 2 , Ri-i) to get the current net- 
work. The algorithm returns a tuple, in which the first component is a boolean that tells 
whether \fi\ = 1; it also returns the new network M(si, s 2 ,Ri), flow /, and distance function 
di to be reused in the next iteration. Smf is inspired by the parametric maximum algorithm 
in [18]. A variant of Smf is used in the first iteration, shown in lines I2.1TT42.131 

This algorithm for sequence of maximum flow problems is called in an improved version 
of SimRel s shown as Algorithm [3j Lines [3^21 13.71 contain the first iteration, very similar to 
the first iteration of Algorithm [1] (lines ll.4Hl.7p . At line 13.41 we prepare for later iterations 
the set 

Listener ( Su s 2 ) = {( u ii u 2) \ u\ € pre(s\) A u 2 G pre(s 2 ) f\ L{u\) = L(u 2 )} , 

where pre(s) = {t G S \ ~P(t,s) > 0}. This set contains all pairs (ui,u 2 ) such that 
the network J\f{u\ , u 2 , R\ ) contains the edge (51,52). Iteration i (for i > 1) of the loop 
(lines [3TT0H3.18P calculates Ri+i from R^. In lines I3.11~H3.141 we collect edges that should 

be removed from N(u\, u 2 , Ri-\) in the sets Dj^ W2 \ At line 13.161 the algorithm Smf 
constructs the maximum flow for parameters using information from iteration % — 1. It uses 
the set D^f 2 ^ to update the network J\f(s\, s 2 , Ri-i), flow a distance function 
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SimRel^ ps (X) 
3.1: Ri <- {(si, s 2 ) E S x S \ L(si) = L(s 2 )} and i <- 1 
3.2: i? 2 «- 

3.3: for all (si, s 2 ) G i?i do 

3.4: Listener ( Sl)S2 ) <— {(iti,^) | Wi G pre(si) A «2 G pre(s 2 ) A L(m) = £(^2)} 

3.5: (mate/i,AA(s 1 ,s 2 ,i?i),/ 1 (si ' S2) ,4 Sl ' S2) ) <~ SMF imf (1, Sl , s 2 , #1) 
3.6: if match then 

3.7: i? 2 ^ J R 2 U{(si,s 2 )} 

while ^ i?j do 
j <- i + 1 

<- and <- \ -Ri 

for all (si,s 2 ) G i?j do 



n (si,s 2 ) 

u %-\ 



for all (si,S2) G fj-i, (^1,^2) G Listener ( Sl)S2 ) H do 



3.8 
3.9 
3.10 
3.11 
3.12 
3.13 

3.14: ^)^^U {(.^..so)} 

3.15 
3.16 

3.17 
3.18 
3.19 



for all (si,S2) G Ri do 

- Smf(z,^( Si , S2 ,^-i) ,fh S2 \d£f 2 \ DtT ] ) 
if match then 

R l+ i <- U {(si,s 2 )} 

return i?. 

Algorithm 3: Improved algorithm for deciding strong simulation for FPSs. 



then it constructs the maximum flow fi for the network M(s\, s 2 , Ri). If Smf returns true, 
(si,s 2 ) is inserted into Ri+\ and survives this iteration (line 13.181) . 

Consider the algorithm Smf and assume that i > 1, At lines 12. 1112.41 we remove the 
edges -Di-i from the network N(s±, s 2 , Ri-i) and generate the preflow fi based on the flow 
fi—l, which is the maximum flow of the network M(s\, s 2 , i?i-i), by 

• setting fi(ui,U2) = for all deleted edges (ui,u 2 ) G Dj-i, and 

• reducing fi(va,\) such that the preflow fi becomes consistent with the (relaxed) flow 
conservation rule. 

The excess e(v) is increased if there exists (v,w) G such that fi^\(v,w) > 0, and 

unchanged otherwise. Hence, fi after line 12.41 is a preflow. The distance function = di 
is still valid for this preflow, since removing the set of edges does not introduce new 

residual edges. This guarantees that, at line 12.51 the preflow algorithm finds a maximum 
flow over the network Af(s\, s 2 , Ri). In line 12.61 Smf returns whether the flow has value 1 
together with information to be reused in the next iteration. (If |/fe| < 1 at some iteration k, 
then \fj\ < 1 for all iterations j > k because deleting edges does not increase the maximum 
flow. In that case, it would be sufficient to return false.) We prove the correctness and 
complexity of the algorithm Smf: 

Lemma 4.3. Let (si,s 2 ) G R±. Then, SMFj n jf returns true iff si s 2 . For some i > 1, 
let M(s±, s 2 , Ri-\), fi-\, and di-\ be as returned by some earlier call to Smf or Smf^. 
Let Di-i = (Ri-i \ Ri) H (post(s\) x post{s 2 )) be the set of edges that will be removed from 
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the network J\f(si, S2, Ri-i) during the (i — \)th call o/Smf. Then, the (i — l)th call o/Smf 
returns true iff si S2- 

Proof. By Lemma |4TT| Smf^ returns true iff |/i| = 1, which is equivalent to si S2- Let 
i > 1. As discussed, at the beginning of line 12,51 the function is a flow (thus a preflow) 
with value 1, and the distance function is a valid distance function. It follows directly 
from the correctness of the preflow algorithm [2] that after line 12.51 fi is a maximum flow 
for M(si, S2, Ri). Thus, Smf returns true (i.e. = 1) which is equivalent to s\ S2- □ 

Lemma 4.4. Consider the pair (s\,S2) and assume that \post(s\)\ < \post(s2)\- -All calls 
to Smf(i, M(sx, S2, •))"' ) related to (s 1; S2) together run in time O{\post{s\)\ \post(s2)\ ). 

Proof. In the bipartite network Af(s\, S2,Ri), the set of vertices are partitioned into subsets 
V\ = pos£(si)U{\} and V2 = post{s2)U{T} as described in Section |4*TT1 Generating the initial 
network (line 12.111) takes time in d?(|Vi| | V2 1 ) . In our sequence of maximum flow problems, 
the number of (nontrivial) iterations, denoted by k, is bounded by the number of edges, i. e., 
k < \E\ < \Vi\ I V2 1 — 1. We split the work being done by all calls to Smf(«, J\f(s%, S2, ■), ■ ■ •) 
together with the initial call to the preflow algorithm (line 12.121 and line 12.50 into edge 
deletions, relabels, non-saturating pushes, saturating pushes. (A non-saturating push along 
an edge (u, v) moves all excess at u to v; by such a push, the number of active nodes never 
increases.) 

All edge deletions take time proportional to Y2i=i l-^«l> w hich is less than the number 
of edges in the network. Therefore, edge deletions take time 0([Vij|V2|). For all v G V, it 
holds that di + i(v) = ck(v), i.e., the labelling function at the beginning of iteration i + 1 is 
the same as the labelling function at the end of iteration i. 

We discuss the time for relabelling and saturating pushes [2]. For a bipartite network, 
the distance of the source can be initialised to d{J~) = 2 \V\ \ instead of |V[, and d(v) never 
grows above 4 | V\ \ for all v E V. For v G V, let I(v) denote the set of nodes containing w such 
that either (v, w) G E or (w, v) G E. Intuitively, it represents edges which could be admis- 
sible leaving v. The time for relabel operations with respect to node v is thus (4 |Vi|)|/(u)|. 
Altogether, this gives the time for all relabel operations: Ylvevd^ l^-l) 1-^(^)1) e ^(l^il l-^D- 
Between two consecutive saturating pushes on (v,w), the distances d(v) and d(w) must in- 
crease by 2. Thus, the number of saturating pushes on edge (v,w) is bounded by 4|Vi|. 
Summing over all edges, the work for saturating pushes is bounded by C?(|Vi| \E\). 

Now we discuss the analysis of the number of non-saturating pushes, which is very 
similar to the proof of Theorem 2.2 in [22] where Max-d version of the algorithm is used. 
Assume that in iteration I < A; of Smf, the last relabelling action occurs. In the Max-d 
version [22], always the active node with the highest label is selected, and once an active 
node is selected, the excess of this node is pushed until it becomes 0. This implies that, 
between any two relabel operations, there are at most n active nodes processed (otherwise 
the algorithm terminates and we get the maximum flow). Also observe that at each time 
an active node is selected, at most one non-saturating push can occur, which implies that 
there are at most n non-saturating pushes between node label increases. Since di{v) is 
bounded by 4|Vi|, the number of relabels altogether is bounded by C(|Vi||V|). Thus, the 
number of non-saturating pushes before the iteration I is bounded by C(|T^||V| 2 ). Since 
the distance function does not change after iteration I any more, inside any of the iterations 
V > Z, there are again at most n — 1 non-saturating pushes. Hence, the number of non- 
saturating pushes is bounded by |li||F| 2 + (fc + 1 - l)(\V\ - 1) G C(|li||y| 2 + k\V\). 
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Since k < | Vi [ | V2 1 — 1, and |V| < 2 1 V2 1 , thus, the overall time complexity amounts to 
C(|Fi||F 2 | 2 ) = O{\post(si)\ \post{s 2 )\ 2 ) as required. □ 

Now we give the correctness and complexity of the algorithm SimRel for FPSs: 

Theorem 4.5. If SimRel fps (.M) terminates, the returned relation equals ^m- 

Proof. By Lemma [4,31 Smf^ (i, s±, S2, Ri) returns true in iteration i = 1 iff si S2; 
Smf(z, M(si, S2, Ri-i), ■ ■ •) returns true in iteration i > 1 iff s% ^r. S2- The rest of the 
correctness proof is the same as the proof of Theorem 14.21 □ 

Theorem 4.6. The algorithm SimRel fps (.M) runs in time 0(m 2 n) and in space 0(m 2 ). 
If the fanout is bounded by a constant, it has complexity 0(n 2 ), both in time and space. 

Proof. We first show the space complexity. In most cases, it is enough to store information 
from the previous iteration until the corresponding structure for the current iteration is 
calculated. The size of the set Listener \ Sl ,s 2 ) ls bounded by |pre(si)| |pre(s2)| where pre(s) = 
{t € S I P(M) > 0}. Summing over all (s 1 ,s 2 ), we get Y^ Sl eS ^ S2 eS |P re ( s i)l |P re ( s 2)| = 
m 2 . Assume we run iteration i. For every pair (si, s 2 ), we generate the set D^l^ and the 
network M(si, s 2 , Ri) together with fi and d{. Obviously, the size of D^J-^ 2 ^ is bounded by 
\post(si)\ \post(s 2 )\. Summing over all (si,s 2 ), we get the bound 0(m 2 ). The number of 
edges of the network M(s\, s 2 , R\) (together with fi and di) is in O{\post(s\)\ \post(s 2 )\). 
Summing over all (s\, s 2 ) yields a memory consumption in 0{m 2 ) again. Hence, the overall 
space complexity is 0{m 2 ). 

Now we show the time complexity. We observe that a pair (s\, s 2 ) belongs to Di in at 
most one iteration. Therefore, the time needed in lines 13.11113. 141 in all iterations together is 
bounded by the size of all sets Listener ( SliS2 ), which is 0(m 2 ). We analyse the time needed 
for all calls to the algorithm Smf. Recall that the fanout g equals max sg s \post(s)\, and 
therefore \post(si)\ < g for i = 1,2. By Lemma |4.4[ the complexity attributed to the pair 
(si, s 2 ) is bounded by 0(g \post{s\)\ \post(s 2 )\). Taking the sum over all possible pairs, we 
get the bound gm 2 € 0(m 2 n). If g is bounded by a constant, we have m < gn, and the 
time complexity is gm 2 < <? 3 n 2 E 0(n 2 ). In this case the space complexity is also 0(n 2 ). □ 

Strong Simulation for Markov Chains. We now consider DTMCs and CTMCs. Since each 
DTMC is a special case of an FPS the algorithm SimRel fps applies directly. 

Let M. = (S, R, L) be a CTMC. Recall that s s' holds if s ~^ e mb(M) s> i n the 
embedded DTMC, and s' is faster than s. We can ensure the additional rate condition by 
incorporating it into the initial relation R. More precisely, initially R contains only those 
pair (s,s') such that L(s) = L(s'), and that the state s' is faster than s, i.e., we replace 
line 13.11 of the algorithm by 

Ri «- {(si, s 2 ) e S x S I L(si) = L(s 2 ) A R(si, S) < R(s 2 , S)} 

to ensure the additional rate condition of Definition ^. 101 In the refinement steps afterwards, 
only the weight function conditions need to be checked with respect to the current relation 
in the embedded DTMC. Thus, we arrive at an algorithm for CTMCs with the same time 
and space complexity as for FPSs. 

Example 4.7. Consider the CTMC in the left part of Figure[7](it has 10 states). Consider 
the pair (si,s 2 ) € R\. The network ftf(s%, S2, Ri) is depicted on the right of the figure. 
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Figure 7: A CTMC example and its network Af(s\, S2,Rx)- 



Assume that we get the maximum flow fx which sends | amount of flow along the path 
T i U2, uZ, \ and ^ amount of flow along /, u\,u^, \. Hence, the check for (s\, S2) is successful 
in the first iteration. The checks for the pairs (ui,us), {u\, U4) and (u2, U3) are also successful 
in the first iteration. However, the check for the pair (1*2,^4) fails, as the probability to go 
from U4 to (73 in the embedded DTMC is |, while the probability to go from 112 to q\ in the 
embedded DTMC is 1. 

In the second iteration, the network J\f(sx, S2, R2) is obtained from Af(si, S2, R\) by 
deleting the edge (u2,u£)- In M(sx, S2, R2), the flows on (u2,Ui) and on \) are set to 
0, and the vertex U2 has a positive excess \. Applying the preflow algorithm, we push the 
excess from U2, along 7Z3, ui,U4 to \. We get a maximum flow /2 for Af(s\, S2, R2) which 
sends \ amount of flow along the path /~,U2,u^,\ and \ amount of flow along f~ ,ux,lu,\- 
Hence, the check for (sx, S2) is also successful in the second iteration. Once the fix-point is 
reached, R still contains (si,^)- 

4.3. Strong Simulation for Probabilistic Automata. In this subsection we present 
algorithms for deciding strong simulations for PAs and CPAs. It takes the skeleton of the 
algorithm for FPSs: it starts with a relation R which is coarser than ^, and then refines R 
until ^ is achieved. In the refinement loop, a pair (s, s') is eliminated from the relation R 
if the corresponding strong simulation conditions are violated with respect to the current 
relation. For PAs, this means that there exists an a-successor distribution \x of s, such that 
for all a-successor distribution // of s', we cannot find a weight function for (/u, //) with 
respect to the current relation R. 

Let A4 = (S, Act, P, L) be a PA. We aim to extend Algorithm [3] to determine the strong 
simulation on PAs. For a pair (si,^), assume that L(si) = L{s2) and that Act(s\) C 
Act(s2), which is guaranteed by the initialisation. We consider line 13.171 which checks the 
condition P(si, •) P(s2 5 •) using Smf. By Definition 12 . 1 1 1 of strong simulation for PAs, 
we should instead check the condition 

Va G Act. Vsi —> fi\. 3s2 — > ^2 with fi± ^2 (4-1) 

Recall the condition /ii \i2 is true iff the maximum flow of the network M{p,Xi ^2, Ri) has 
value one. Sometimes, we write M(sx, a, fJ>x, S2, fi2,Ri) to denote the network N(fix, Ri) 
associated with the pair (si, S2) with respect to action a. 

Our first goal is to extend Smf to check Condition 14.11 for a fixed action a and a- 
successor distribution fix of s%. To this end, we introduce a list Sim^ Sl,a,tJ,1,S2 ^ that contains 
all potential candidates of a-successor distributions of S2 which could be used to establish 
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ActSmf^, Sirrii-uAfijiuHz, Ri-i), fi-i, di-i, A-i) 

4.1: Simi <— Sim. L -i 

4.2: (match, M(m, /i 2 , -Ri), /i, d») <- Smf(«, AT(/ii, /i 2 , i^-i), /i_i, A-i) 

4.3: if match then 

4.4: return (true, Simi, pL 2 , Rt), ft, (k) 

4.5: Simi <— tail(5irrij) 

4.6: while -lempty (Simi) do 

4.7: ^2 head(S'imj) 

4.8: (match,J\f(m,iJ,2,Ri),fi,di) <- SMF init (i, fj,t, /j, 2 , R4) 
4.9: if match then 

4.10: return (true, Sirrii,N '(/ii, //2j -Ri), fi,d{) 

4.11: Szmj <— tail(S'imi) 

4.12: return (false,®, NIL, NIL, NIL) 

ACTSMF«(i,/ii, Simi,Ri) 
goto line 14.61 

Algorithm 4: Subroutine to calculate whether si s 2) as far as si — > /xi is concerned. The 
parameter Sim denotes the subsets of a-successor distributions of s 2 serving as candidates 
for possible \i 2 . 

the condition \i\ C fi fj, 2 for the relation R considered. The set Sim^ 1 ' 01 '^ 1 ' 82 ' is represented 
as a list. This and some subsequent notations are similar to those used by Baier et al. in [3]. 
We use the function head(-) to read the first element of a list; tail(-) to read all but the first 
element of a list; and empty (•) to check whether a list is empty. As long as the network 
for a fixed candidate fi 2 = head(S , im^ 1,a '' 11 ' S2 ^) allows a flow of value 1 over the iterations, 
we stick to it, and we can reuse the flow and distance function from previous iterations. If 
by deleting some edges from M(pL%, \i 2 , R), its flow value falls below 1, we delete \i 2 from 
Sj m (si,a,[ii,s 2 ) an( j nex t candidate. 

The algorithm ActSmf, shown as Algorithm HJ implements this. It has to be called 
for each pair (si,s 2 ) and each successor distribution s\ — > fii of s±. It takes as input the 
list of remaining candidates Simf^f i, ^ 1,S2 \ the information from the previous iteration (the 
network N(m, fJ, 2 , Ri-i), flow fi—i, and distance function dj-i), and the set of edges that 
have to be deleted from the old network A-i- 

Lemma 4.8. Let (s\,s 2 ) £ R\, a E Act(s\), and [i\ such that s\ — > /ii. Let Sim\ = 
Steps a (s 2 ). Then ActSmFj^ returns true iff 3[i 2 with s 2 — * \i 2 A /Ui jjh 2 . For some 
i > 1, let Simi— i, M(n\, \i 2 , Ri-i), and be as returned by some earlier call to 
ActSmf or ActSmf^. Let Di~\ = (R4-1 \ R4) n (Supp(^i) x Supp(/j, 2 )) be the set of 
edges that will be removed from the network during the (i — l)th call of ActSmf. Then, the 
(i — l)th call of algorithm ActSmf returns true iff: 3fi 2 with s 2 — > fi 2 A [i\ Qfy fi 2 . 

Proof. Once Smf returns false because the maximum flow for the current candidate [i 2 has 
value < 1 , it will never become a candidate again, as edge deletions cannot lead to increased 
flow. The correctness proof is then the same as the proof of Lemma 14.31 □ 
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SimRel fa (M) 



5.1: R\ <- {(si, s 2 ) e S x S \ L(si) = L(s 2 ) A Act(si) C ^(s^} and i <- 1 

5.2: R 2 «- 

5.3: for all (si,,s 2 ) £ -Ri do 

5.4: Izstener( SljS2 ) <— {(it a , a, fix, u 2 , fi 2 ) \ L(ux) = L(u 2 ) f\ux ^ ftiAn 2 fi 2 

A/ii(si) > A/X 2 (s 2 ) > 0} 
5.5: for all a G -Aci(si), /xi G Steps a (sx) do 

5.6: Sim^ u0 ^ u82) «- Sieps a ( S2 ) 

5.7: (raafcfc^ , Szm^ 1 ' W2 \ A/"(si, a, Ml , s 2 , /x 2 , #1), /f c , df c ) 

<- AcTSMF mii ( 1 , (xx , ^ml Sl ' S2) ,Rx) 

5 - 8: if A ae Arf( Sl ) A we step Sa ( sl ) match^ then 

5.9: i? 2 <- R 2 U {(si,s 2 )} 

5.10: while R l+1 / i?j do 
5.11: + l 

5.12: i? m <- and A-i <- -Rj-i \ Ri 

5.13: for all (si,s 2 ) G i?, a € Aci(si), /xi G Steps a (sx), fi 2 G Steps a {s 2 ) do 
5.14: £ ,(«wx,»2.A«) ^ 

5.15: for all (si,s 2 ) G Di-x, {ux,a, fix,u 2 , fi 2 ) G Listener ( SljS2 ) ^° 
5.16: if (ui,u 2 ) G -Ri-i then 

5.17: ^(ui.a.w^.w) ^ £) («x,a,/i 1 ,«a,/i a ) (j {( Sl> S2 )} 

5.18: for all (si, s 2 ) G i? do 

5.19: for all a G j4ct(si), /xi G Steps a (sx) do 

5.20: (match^, Sim[ ai ' a ^ u82 \ Af(sx, a, fix,s 2 ,fi 2 , Ri), /f c , df c ) 

<- ActSmf(z, Sim^{ a ' m > S2 \ M{ Sl ,a,fj,x,s 2 ,Li 2 ,Ri-x), 

fare mrc rtarc \ 

5 - 21: if f\aeAct( sl ) A w eSfcjw«(»i) ™atch a ^ then 

5.22: Ri+X^ Ri +1 U{(sx,S 2 )} 

5.23: return Ri 



Algorithm 5: Algorithm for deciding strong simulation for PAs, where arc denotes the 
associated parameter (sj, a, /xi, s 2 , fi 2 ). 

The algorithm SimRel pa for deciding strong simulation for PAs is presented as Algo- 
rithm [5j During the initialisation (lines ETHESl intermixed with iteration 1 in lines 15. 7l45.9p . 
for (si,s 2 ) G i?i and sx —> fix, the list Sim^ 1 ' 01 '^ 1 is initialised to Steps Q (s 2 ) (line 15 .6[) . 
as no a-successor distribution can be excluded as a candidate a priori. As in SimRel fps , 
the set Listener ( Sl)S2 ) for (sx,s 2 ) is introduced which contains tuples (u\, a, fix, u 2 , /x 2 ) such 
that the network M(ux,a, fix,u 2 , fi 2 , R\) contains the edge (si,^). 

The main iteration starts with generating the sets £)^^ a >^ 1 ' U2ilJ ' 2 ' [ n lines I5.13h5.17I in 
a similar way as SimRel fps . Lines check Condition 14.11 by calling ActSmf for 

each action a and each a-successor distribution fix of sx- The condition is true if and only 
if matcha^ is true for all a G Act(sx) and fix G Steps a (sx)- In this case we insert the pair 
(si,s 2 ) into (line 15.221) . We give the correctness of the algorithm: 

Theorem 4.9. When SimRel pa (A4) terminates, the returned relation equals ~^m- 
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Proof. The proof follows the same lines as the proof of the correctness of SimRel s in 
Theorem 14.51 The only new element is that we now have to quantify over the actions and 
successor distributions as prescribed by Definition 12.111 This translates to a conjunction in 
lines 15.81 and I5.2T1 of the algorithm. Exploiting Lemma 14.81 we get the correctness. □ 

Now we give the complexity of the algorithm: 

Theorem 4.10. The algorithm SimRel^ a (A1) runs in time 0(m 2 n) and in space 0{m 2 ). 
If the fanout of M. is bounded by a constant, it has complexity 0{n 2 ), both in time and 
space. 

Proof. We first consider the space complexity. We save the sets the networks 

J\f(si,a, fix, S2, fi2, R%) which are updated in every iteration, Listener ( sliS2 ) and the sets 

Sirnf 1,a ' IJ ' l,s '^ . The size of the set jj^ 1 ' 01 '^ 1 ' 82 '^ j s \ n 0(\fii\ \fi2\), which is the maximal 
number of edges of M(sx, a, fix, S2,fi2,Ri)- Summing over all (sx, a, fix, s%, ^2), we get: 

EE E E E \^\<™ 2 (4.2) 

siGS a&Act{s\) fii£Steps a (si) S2GS ii2&Steps a (s2) 

Similarly, the memory needed for saving the networks has the same bound 0(m 2 ). Now 
we consider the set Listener ( SljS2 ) for the pair (si,S2) £ R\. Let (u±, a, fix, U2, 1^2) £ 
Listener ( Sl ,s 2 )- Then, it holds that s\ E Supp(fix) and S2 € Supp(fi2). Hence, the tuple 
(iti, a, fix, U2, H2) can be an element of Listener \ Sl ,s 2 ) °f some arbitrary pair (s±, S2) at most 
I Mi I I 1 times, which corresponds to the maximal number of edges between the set of nodes 
Supp(ni) and Supp(fi2) m N(s\, a, fix, S2, fJ-2, Ri)- Summing over all (s\,a, fi\,S2, ^2)-, we 
get that memory needed for the set Listener is also bounded by 0(m 2 ). For each pair (si, S2) 

and si ^ fix, the set Sim^ 1 '"'^ 1 ' 82 ^ has size \Steps a {s2)\- Summing up, this is smaller than 
or equal to m 2 according to Inequality 14.21 Hence, the overall space complexity amounts to 
0(m 2 ). 

Now we consider the time complexity. All initialisations (lines [5TTH5.6l of SimRel^ a and 
the initialisations in AcTSMFj n j(, which calls SMFj ra jt) take 0(m 2 ) time. We observe that a 
pair (sx,S2) belongs to Di during at most one iteration. Because of the Inequality 14.21 the 
time needed in lines 15.13115. 171 is in 0(m 2 ). The rest of the algorithm is dominated by the 
time needed for calling Smf in line 14.21 of ActSmf. By Lemma 14.41 the time complexity 
for successful and unsuccessful checks concerning the tuple (si, a, fix, S2, ^2) is bounded by 
0(g\fix\ I/X2I). Taking the sum over all possible tuples (si, a, fix, $2, ^2) we get the bound 
gm 2 according to Inequality 14.21 Hence, the complexity is 0(m 2 n). If the fanout g is 
bounded by a constant, we have m < gn. Thus, the time complexity is in the order of 
0(n 2 ). In this case the space complexity is also 0(n 2 ). □ 

Remark 4.11. Let M = (S,Act,P,L) be a PA, and let k = J2 s esY^ a <EAct(s) \Steps a (s)\, 
called the number of transitions in [3], denote the number of all distributions in M. The 
algorithm for deciding strong simulation introduced by Baier et al. has time complexity 
0((kn 6 + k 2 n 3 )/ log n), and space complexity 0(k 2 ). The number of distributions k and 
the size of transitions m are related by k < m < nk. The left equality is established if \fi\ = 1 
for all distributions, and the right equality is established if \fi\ = n for all distributions. 

The decision algorithm for strong simulation for CPAs can be adapted from SimRel^ a 
in Algorithm [5] easily: Notations are extended with respect to rate functions instead of 
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distributions in an obvious way. To guarantee the additional rate condition, we rule out 
successor rate functions of S2 that violate it by replacing line 15.61 by: 

For each pair (si,S2)> and successor rate functions rj G Steps a (si) (i = 1,2), the subrou- 
tine for checking whether r\ C/j. T2 is then performed in the network A/"(/u(ri), [i{r<i), Ri). 
Obviously, the so obtained algorithm for CPAs has the same complexity 0(m 2 n). 

4.4. Strong Probabilistic Simulation. The problem of deciding strong probabilistic sim- 
ulation for PAs has not been tackled yet. We show that it can be computed by solving LP 
problems which are decidable in polynomial time [27J. In Subsection 14.4. 1] , we first present 
an algorithm for PAs. We extend the algorithm to deal with CPAs in Subsection 14.4.21 

4.4.1. Probabilistic Automata. Recall that strong probabilistic simulation is a relaxation 
of strong simulation in the sense that it allows combined transitions, which are convex 
combinations of multiple distributions belonging to equally labelled transitions. Again, the 
most important part is to check whether s\ ^ R S2 where R is the current relation. By 
Definition 12.151 it suffices to check L(si) = L(s2) and the condition: 

Va G Act. Vsi ^* n\. 3&2 ~» ^2 with /ii ^2 (4-3) 

However, since the combined transition involves the quantification of the constants c, G 
[0, 1], there are possibly infinitely many such fj,2- Thus, one cannot check fi\ Cr fi2 for each 
possible candidate fj,2- The following lemma shows that this condition can be checked by 
solving LP problems which are decidable in polynomial time [27} [5B] . 

Lemma 4.12. Let M = (S,Act,P,L) be a given PA, and let RQ S x S. Let (si, S2) 6 R 
with L(s\) = L{s2) and Act{s\) C Act{s2). Then, s\ ^ R S2 iff for each transition s\ —> fi, 
the following LP has a feasible solution: 

k 

i=l 

0<Ci<l V t = 1, . . . , Jfe (4.5) 
0</ (s , t) <l V(a,t)€R± (4.6) 

M*) = E A-,*) ^seS ± (4.7) 
teR x (s) 

k 

E /(-,*) =E Ci ^(*) VteS ± (4.8) 
where k = \Steps a (s2)\ > and Steps a {s2) = {/ii, . . . , fJ>k}- 

Proof. First assume that s\ ^ R 82- Let s\ — > By the definition of simulation up to 
R for strong probabilistic simulation, there exists a combined transition S2 ~> ^ c with 
M E-R. Mc- Let Steps a {s2) = {Ml> • • • i Mfc} where = | Steps a ($2) | • Now Act(si) C Aci(s2) 
implies fc > 0. By definition of combined transition (Definition I2.14|) . there exist constants 
ci, . . . , Ck G [0, 1] with Yli=i °i = 1 such that fj, c = Ym=i CitH- Thus Constraints 14.41 and 14.51 
hold. Since fi fi c , there exists a weight function A : S± x S± — > [0, 1] for /x c ) with 
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SimRel™' p (7W) 

6.1: Ri <- {(si, s 2 ) e S x S \ L(si) = L(s 2 ) A Act(si) C ,4ci(s 2 )} and j <- 

6.2: repeat 

6.3: i <— « + 1 

6.4: i? m <- 

6.5: for all (si,s 2 ) € i?j do 

6.6: for all a G Act{s\), [ii G SYeps^si) do 

6.7: matcha^i <— LP(si,a, fii, S2) 

6 - 8: if AaeAct( Sl ) A w eSte P 5 (5 1 ) match atlil then 

6.9: i? m <- i^+i U {(s!,s 2 )} 

6.10: until = i?j 

6.11: return i?j 

Algorithm 6: Algorithm for deciding strong probabilistic simulation for PAs. 

respect to R. For every pair (s,t) G R±, let fr 8t t) := A(s,i). Thus, Constraint 14.61 holds 
trivially. By Definition 12.71 of weight functions, it holds that (i) A(s,t) > implies that 
(s,t) G R±, (ii) n(s) = J2tes ± A ( s >0 for s e s ±, and ( ni ) vS) = J2 s eS ± for a11 

t G Sj_. Observe that (i) implies that for all (s,i) G" we have that A(s,t) = 0. Thus, 
(ii) and (hi) imply Equations 14.71 and 14.81 respectively. 

Now we show the other direction. Let k = \Steps a (s2)\ and Steps a (s2) = {Ml> • • • 
By assumption, for each s\ —> fi, we have a feasible solution ci, . . . , c& and fr s>t \ for all (s, t) G 

which satisfies all of the constraints. We define \i c = Y^i=i c i^i- By Definition 12.141 
[i c is a combined transition, thus s 2 ~* Mc- It remains to show that /j, /i c . We define a 
function A as follows: A(s, t) equals /( s t ) if (s,t) G R± and otherwise. With the help of 
Constraints 14.61 14-71 and 14.81 we have that A is a weight function for (/x, // c ) with respect to 
R, thus fj, Qr fi c . □ 

Now we are able to check Condition 14.31 by solving LP problems. For a PA A4 = 
(S, Act,P, L), and a relation R C S X S, let (si,s 2 ) G i? with L(si) = £(s 2 ) and Aci(si) C 
Arf(s 2 ). For s± —> Hi, we introduce a predicate LP(si, a, fj,, s 2 ) which is true iff the LP 
problem described as in Lemma 14.121 has a solution. Then, si 7^ s 2 iff the conjunc- 
tion A Q eArf( S i) h^esteps^s!) LP ( s ii a > Ml) ^2) is true. The algorithm, which is denoted by 
SiMREL^ A>p (A / i), is depicted in Algorithm El It takes the skeleton of SimRel s (.M). The key 
difference is that we incorporate the predicate LP{s\, a, /xi, s 2 ) in line 16.71 The correctness 
of the algorithm SimRel^ A ' p (.M) can be obtained from the one of SimRel s (A1) together 
with Lemma 14.121 We discuss briefly the complexity The number of variables in the LP 
problem in Lemma f4.12l is /e+|i?|, and the number of constraints is l+fc+|i?|+2 \S\ G C(|i?|). 
In iteration i of SimReLs A ' p (.M), for (s\, s 2 ) G Ri and si —> fJ-i, the corresponding LP prob- 
lem is queried once. The number of iterations is in 0(n 2 ). Therefore, in the worst case, one 
has to solve n 2 J2ses J2aeAct(s) £ M eSk?w«,(«) 1 e 0(n 2 m) many such LP problems and each 
of them has at most 0(n 2 ) constraints. 

4.4.2. Continuous-time Probabilistic Automata. Now we discuss how to extend the algo- 
rithm to handle CPAs. Let A4 = (S, Act, R, L) be a CPA. Similar to PAs, the most impor- 
tant part is to check the condition s\ r^p s 2 for some relation R C S x S. By Definition 12. 191 



28 



L. ZHANG, H. HERMANNS, F. EISENBRAND, AND D. N. JANSEN 



SimRei^ 

7.1: i?i <- {(si, s 2 ) eSxS L(si) = L(s 2 ) A Act(si) C ,4ci(s 2 )} and i <- 

7.2: repeat 

7.3: i <— i + 1 

7.4: i? m <- 

7.5: for all (si,s 2 ) G ^ d° 

7.6: for all a G Act(s\), r\ G Steps a (s{), E £ E(s 2 ) do 

7.7: match ajri) E LP'(si,a, r±, s 2 ,E) 

7 - 8: if AaeAci( Sl ) Ane^ep^fsi) /\eee{s 2 ) match a} r ltE then 

7.9: ^-Rj+iU{(s 1 ,s 2 )} 

7.10: until R i+ i = Ri 

7.11: return R{ 

Algorithm 7: Algorithm for deciding strong probabilistic simulation for CPAs. 

it suffices to check L(si) = L(s 2 ) and the condition: 

Va G Act. Vsi — » ri. 3s 2 ~> r 2 with /i(ri) C fi /u(r 2 ) A ri(S') < r 2 (5) (4-9) 

Recall that for CPAs only successor rate functions with the same exit rate can be combined 
together. For state s G S, we let E(s) := {r(S) \ s —> r} denote the set of all possible exit 
rates of a-successor rate functions of s. For E G E(s) and a G Act(s), we let Steps^(s) = 
{r G Steps a {s) \ r(S) = E} denote the set of a-successor rate functions of s with the same 
exit rate E. As for PAs, to check the condition s% s 2 we resort to a reduction to LP 
problems. 

Lemma 4.13. Let M = (S, Act, R, L) be a given CPA, and let R C S x S. Let (s\, s 2 ) G R 
with L(si) = L(s 2 ) and that Act{s\) C Act(s 2 ). Then, s\ s 2 iff for each transition 

s± — > r either r(S) = 0, or there exists E G E(s 2 ) with E > r(S) such that the following 
LP has a feasible solution, which consists of Constraints \4-4\ \4-5\ \4-6\ °f Lemma \4-l%\ and 
additionally: 

r(s)=r(S) f(s,t) VsGS ± (4.10) 
teR±(s) 

k 

E E /(m) = E c ^w vtG ^ ( 4 - n ) 

where k = \ Steps^ (s)| with Steps ^ (s) = {ri, . . . , r^}. 

Proof. The proof follows the same strategy as the proof of Lemma f4.12| in which the induced 
distribution of the corresponding rate function should be used. □ 

Now we are able to check Condition 14.91 by solving LP problems. For a CPA A4 = 
(S, Act, II, L), and a relation R C S x S, let (si,s 2 ) G R with L(si) = L(s 2 ) and Act(s\) C 
Aci(s 2 ). For si n, and E G E(s 2 ), we introduce the predicate LP'{s\, a, r\, s 2 , E) which 
is true iff E > r\(S) and the corresponding LP problem has a solution. Then, s\ s 2 iff 
the conjunction f\ aeAct(si) K ieS teps a (s 2 ) AfieS^) LP'{s u a, n, s 2 , E) is true. The decision 
algorithm is given in Algorithm As complexity we have to solve 0(n 2 m) LP problems 
and each of them has at most 0(n 2 ) constraints. 
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5. Algorithms for Deciding Weak Simulations 

We now turn our attention to weak simulations. Similar to strong simulations, the core of 
the algorithm is to check whether s± ^.r S2, i.e., S2 weakly simulates s% up to the current 
relation R. As for strong simulation up to R, s± ^.r S2 does not imply s\ S2, since 
no conditions are imposed on pairs in R different from (s\,S2)- By the definition of weak 
simulation, for fixed characteristic functions 5i (i = 1,2), the weight function conditions 
can be checked by applying maximum flow algorithms. Unfortunately, ^-functions are not 
known a priori. Inspired by the parametric maximum flow algorithm, in this chapter, we 
show that one can determine whether such characteristic functions <5j exist with the help of 
breakpoints, which can be computed by analysing a parametric network constructed out of 
P(si, -),P(s2, •) and R. We present dedicated algorithms for DTMCs in Subsection 15.11 and 
CTMCs in Subsection EZS 

5.1. An Algorithm for DTMCs. Let M = (S,P,L) be a DTMC. Let R C S x S be 

a relation and s\ R S2- Whether S2 weakly simulates s\ up to R is equivalent to whether 
there exist functions 5i : S — > [0, 1] such that the conditions in Definition 12.211 are satisfied. 
Assume that we are given the [/^-characterising functions <5j. In this case, s\ ^Sr S2 can be 
checked as follows: 

• Concerning Condition we check whether for all v £ S with 6\(v) < 1 it holds that 
v R S2- Similarly, for Condition [lb, we check whether for all v 6 S with S^iv) < 1 it 
holds that s\ R v. 

• The reachability condition can be checked by using standard graph algorithms. In more 
detail, for each u with 6\(u) > 0, the condition holds if a state in R(u) is reachable from 
S2 via R(si) states. 

• Finally consider Condition [2j From the given 6{ functions we can compute Ki . In case 
of that K\ > and K2 > 0, we need to check whether there exists a weight function 
for the conditional distributions P(si, -)/K\ and P(s2> ')/-^2 with respect to the current 
relation R. From Lemma |4.1( this is equivalent to check whether the maximum flow for 
the network constructed from (P(si, •)/iTi,P(s2, O/^) an d R has value 1. 

To check s\ ^.r S2, we want to check whether such <5j functions exist. The difficulty is 
that there exist uncountably many possible 5i functions. In this section, we first show that 
whether such Si exists can be characterised by analysing a parametric network in Subsec- 
tion [5]TTJ Then, in Subsection 15.1.21 we recall the notion of breakpoints, and show that the 
breakpoints play a central role in the parametric networks considered: only these points need 
to be considered. Based on this, we present the algorithm for DTMCs in Subsection 15.1.31 
An improvement of the algorithm for certain cases is reported in Subsection 15.1.41 

5.1.1. The Parametric Network J\f(j). Let 7 <G M> . Recall that jV(P(si, •), 7P(s2, •)> R ) 
is obtained from the network A/"(P(si, •), P(s2, •), R) be setting the capacities to the sink \ 
by: c(t,\) = 7P(s2,i). If si,S2,R are clear from the context, we use M(-y) to denote the 
network J\f(P(s\, -),j~P(s2, -),R) for arbitrary 7 G M>o- 

We introduce some notations. We focus on a particular pair (si, S2) € R, where R is the 
current relation. We partition the set post(si) into MUi (for: must be in Ui) and PVi (for: 
potentially in Vi). The set PV\ consists of those successors of s\ which can be either put 
into U\ or V\ or both: PV\ = post{s\) {~\R~ l {s2). The set MU\ equals post{s\)\PV\, which 
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Figure 8: A simple DTMC. 



consists of the successor states which can only be placed in U±. The sets PV2 and MU 2 
are defined similarly by: PV 2 = post(s 2 ) Pi R(si) and MU2 = post(s 2 )\PV 2 . Obviously, 
5i(u) = 1 for u G MUi for i = 1,2. 

Example 5.1. Consider the DTMC in Figure [8] and the relation R = {(s\, S2), {s±, 1)%), 
(vi,s 2 ), (ui,u 2 ), (oi,o 2 ), (01, v 2 ), (fi,o 3 ), (vi,v 2 ), (o 2 , 01)}. We have PV\ = {vi} and PV 2 = 
{v 2 }- Thus, MUi = {ni,oi}, MU 2 = {^2,02,03}. 

We say a flow function / of Affa) is valid for ^(7) iff / saturates all edges {f,u\) with 
u\ G MU\ and all edges (H2, \) with 112 G MU 2 . If there exists a valid flow / for Affa), we 
say that 7 is valid for Affa). The following lemma considers the case in which both si and 
S2 have visible steps: 

Lemma 5.2. Let s± R s 2 . Assume that there exists a state s± G post(si) such that s[ G" 
R^ 1 (s 2 ), and s' 2 G post(s2) such that s' 2 G" R(s±). Then, si ~r S2 iff there exists a valid 7 
for Af{y). 

Proof. By assumption, we have that s[ G MUi for i = 1, 2, thus MC/j 7^ 0, and it holds that 
6i(sfi = 1 for i = 1,2. 

We first show the on/y i/ direction. Assume s\ tSr s 2 , and let Si, Ui, Vi, Ki, A (for 
i = 1,2) as described in Definition 12.211 Since MUi 7^ for i = 1,2, both K\ and -K2 are 
greater than 0. We let 7 = K\/K 2 . For s, s' G 5, we define the function / for Af(^y): 

f(r,s) = P(s 1 ,s)8 1 (s), f(s,t)=K 1 A(s,t), f(s,\)= 1 P( S 2,s)5 2 (s) 

Since 5 l (s) < 1 (i = 1, 2) for s G 5, /(/, s) < P(si, s) and /(*, \) < 7P(s 2 , a). Therefore, / 
satisfies the capacity constraints. / also satisfies the conservation rule: 

/( s ,5) = ifiA( s ,s) = p( Sl)S )5i( s ) =/(;,«) 

=^iA(5,s) =7^ 2 A(5, S )= 7 P(a2,a)<52W = /(»A) 
Hence, / is a flow function for Af(-y). For u\ G MU\, we have <5i(ui) = 1, therefore, 
/(/,ui) = P(si,ui). Analogously, f(u2,\) = "fP(s 2 ,U2) for U2 G MU 2 . Hence, / is valid 
for A/"(7), implying that 7 is valid for Affa). 

Now we show the i/ direction. Assume that there exists 7 > and a valid flow / 
for Af(j). The function <5i is defined by: 6i(s) equals /(/, s)/P(si, s) if s G post(si) and 
otherwise. The function 5 2 is defined similarly: ^(s) equals f(s, \)/jP(s 2 , s) if s G post (32) 
and otherwise. Let the sets C/j and be defined as required by Definition 12.211 It follows 
that 

Ki= £*i(s)p(*i,*) = 52f(r,s) = nr,ui) 

K 2 = 2^ S 2 (S)P(S2, s) = ^ — — = 

sec/2 set/ 2 
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Figure 9: Left: The network M{2) of the DTMC in Figure [HJ Right: The transformed 
network M (2) for Af(2). 



Since the amount of flow out of f is the same as the amount of flow into \, we have 
Ki/K 2 = 7. Since / M{/ 8 C Ui for i = 1, 2, both of Ki and K 2 are greater than 0. We 
show that the Conditions QJi and[Tb of Definition 12.211 are satisfied. For v\ G V%, we have 
that 5\{vi) < 1 which implies that /(/, vi) < P(si, v±). Since / is valid for Af("f), and since 
the edge is not saturated by /, it must hold that v\ G PV\. Therefore, v\ R s 2 . 

Similarly, we can prove that s\ R v 2 for v 2 G V 2 . 

We define A(w,w') = f(w,w')/Ki for w,w' G S. Assume that A(w,w') > 0. Then, 
f(w,w') > 0, which implies that (w,w') is an edge of 7V(7), therefore, (w,w') € R. By the 
flow conservation rule, f(f,w) > f(w,w') > 0, implying that 5i(io) > 0. By the definition 
of Ui, we obtain that w £ U\. Similarly, we can show that w' G U 2 . Hence, the Condition [2al 
is satisfied. To prove Condition | 

AKZ7 a ) 



^ /(u>,u 2 ) _ f(w,U 2 ) (*) /(/,«;) _ Ji(w)P(ai,«;) 



^ r Kl ^! ^ K x 

U2&U2 

where equality (*) follows from the flow conservation rule. Therefore, for w G S we have 
that K\ A(w, U 2 ) = P(si,w)5i(w). Similarly, we can show K 2 A(Ui,w) = P(s 2 ,w)5 2 (w). 
Condition [2b] is also satisfied. As K\ > and K 2 > 0, the reachability condition holds 
trivially, hence, s± s 2 . □ 

Example 5.3. Consider again Example 15.11 with the relation R = {(s±, s 2 ), (si, v 2 ), (v\, s 2 ), 
(ui,u 2 ), (01,02), (01,^2), (yi, 03), (vi,v 2 ), (02, 01)}. The network M(2) is depicted on the left 
part of Figure ES Edges without numbers have capacity 00. It is easy to see that 2 is 
valid for jV(2): the corresponding flow sends \ amount of flow along the path /, u\,U2, \, I 
amount of flow along the path /, 01,02, \, \ amount of flow along the path /',ox,v 2 ,\, and 
I amount of flow along the path v\ , 03 , \, . 

For a fixed 7 G M>o, we now address the problem of checking whether there exists a 
valid flow / for M (7). This is a feasible flow problem (/ has to saturate edges to MU\ 
and from MU 2 ). As we have discussed in Section [3J it can be solved by applying a simple 
transformation to the graph (in time 0{\MU\\ + \MU 2 \)), solving the maximum flow problem 
for the transformed graph, and checking whether the flow saturates all edges from the new 
source. 

Example 5.4. Consider the network M{2) on the left part of Figure GO Applying the 
transformation for the feasible flow problem described in Section [31 we get the transformed 
network A/"t(2) depicted on the right part of Figure EJ It is easy to see that the maximum 
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flow h for A/i (2) has value Namely: It sends | amount of flow along the path f ', ui,u^, \' , 
| amount of flow along the path 01,02, \', | amount of flow along o\,V2, \, A V', | 
amount of flow along \, /, v\, 03, \' , and ^ amount of flow along /',\,/,V- Thus, it uses 
all capacities of edges from /'. This implies that 2 is valid for the network A"(2). 



5.1.2. Breakpoints. Consider the pair (si, S2) € -R- Assume the conditions of Lemma [52] are 
satisfied, thus, to check whether si s 2 it is equivalent to check whether a valid 7 for Af(^) 
exits. We show that only a finite possible 7, called breakpoints, need to be considered. The 
breakpoints can be computed using a variant of the parametric maximum flow algorithm. 
Then, si s 2 if and only if for some breakpoint it holds that the maximum flow for the 
corresponding transformed network At (7) has a large enough value. 

Let \V\ denote the number of vertices of A"(7). Let ^(7) denote the minimum cut 
capacity function of the parameter 7, which is the capacity of a minimum cut of A"(7) as a 
function of 7. The capacity of a minimum cut equals the value of a maximum flow. If the 
edge capacities in the network are linear functions of 7, ^(7) is a piecewise-linear concave 
function with at most [V| — 2 breakpoints [18j . i. e., points where the slope |~ changes. The 
\V\ — 1 or fewer line segments forming the graph of ^(7) correspond to \V\ — 1 or fewer 
distinct minimal cuts. The minimum cut can be chosen as the same on a single linear piece 
of ^(7), and at breakpoints certain edges become saturated or unsaturated. The capacity 
of a minimum cut for some 7* gives an equation that contributes a line segment to the 
function ^(7) at 7 = 7*. Moreover, this line segment connects the two points (71,^(71)) 
and (72,^(72)), where 71,72 are the nearest breakpoints to the left and right, respectively. 

Example 5.5. Consider the DTMC in Figure [HJ together with the relation R = {(si, S2), 
(si,v 2 ), (vi,s 2 ), {ui,u 2 ), (oi,o 2 ), (oi,v 2 ), (vi,o 3 ), (vi,v 2 ), (o 2 ,oi)}. The network M(j) for 
the pair (si,s 2 ) is depicted in Figure [TU1 

There are two breakpoints, namely I and 2. For 7 < I, all edges leading to the sink can 
be saturated. This can be established by the following flow function /: sending ^ amount 
of flow along the path /,«i,U2, \, ^ amount of flow along the path /',o\,o 2 ,\, amount 
of flow along the path /',o\,v 2 ,\, amount of flow over 1/1,03, \, amount of flow 
along the path /',vi,v 2 ,\. The amount of flow out of node 01, denoted by f(o±,S), is j^. 
Given that 7 < I, we have that f(oi,S) < \. Similarly, consider the amount of flow out of 
node vx, which is denoted by f(v\, S), is ^ which implies that f(v\,S) < \. The maximum 
flow thus has value |/| = ^ + ^ + "24" + lis + Ts" = ^' Thus the value of the maximum flow, 
or equivalently the value of the minimum cut, for 7 < | is ^(7) = 7. 
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Figure 11: The value of the maximum flow, or equivalently the value of the minimum cut, 
as a function of 7 for the network in Figure [TUl 

Observe that for 7 = |, the edges to v\ and o\ are saturated, i.e., we have used full 
capacities of the edge (A^i) an d (/,oi). Thus, by a greater value of 7, although the 
capacities c({v2, 02,03}, \) increase (become greater than |), no additional flow can be sent 
through {vi, oi}. For the other breakpoint 2, we observe that for a value of 7 < 2, we can 
still send ^ through the path /, u\, u 2 , \, but for 7 > 2, the edge to u\ keeps saturated, thus 
the amount of flow sent through this path does not increase any more. Thus, for 76 [|, 2], 
the maximum value, or the value of the minimum cut, is | + ^. The first term | corresponds 
to the amount of flow through v\ and o\. The breakpoint | is not valid since the edge to 
u\ can not be saturated. As discussed in Example 15.41 the breakpoint 2 is valid. The curve 
for ^(7) is depicted in Figure [TT1 

In the following lemma we show that if there is any valid 7, then at least one breakpoint 
is valid. 

Lemma 5.6. Assume 7* G (71,72) where 71,72 are two subsequent breakpoints 0/^(7), or 
71 = and 72 is the first breakpoint, or 71 is the last breakpoint and 72 = 00. Assume 7* 
is valid for N{^*), then, 7 is valid for N(-y) for all 7 G [71,72]- 

Proof. Consider the network M (7*). Assume that the maximum flow fy is a valid maximum 
flow for A r (7*). 

Assume first 7' G (7*, 72]. We use the augmenting path algorithm [lj to obtain a 
maximum flow /* in the residual network .A// ,(7'), requiring that the augmenting path 
contains no cycles, which is a harmless restriction. Then, fy := / 7 * + /* is a maximum 
flow in A/"(7 ; ). Since / 7 * saturates edges from / to MU\, fy saturates edges from / to 
MU\ as well , as flow along an augmenting path without cycles does not un-saturate edges 
to MU\. We choose the minimum cut (X, X') for Af(^*) with respect to fy such that 
MU2 H X' = 0, or equivalently Mt7 2 C X. This is possible since / 7 * saturates all edges 
(u2,\) with U2 G MU%. The minimum cut for fy, then, can also be chosen as (X, X'), as 
(7', «(Y)) lies on the same line segment as (7*,k(7*)). Hence, fy saturates the edges from 
MU2 to \, which indicates that fy is valid for J\f{^'). Therefore, 7' is valid for M^') for 
i € (7*,72j- 

Now let 7' G [71, 7*). For the valid maximum flow fy we select the minimal cut (X, X') 
for J\f(j*) such that MU\ D X = 0. Let d denote a valid distance function corresponding 
to fy*. We replace fy(v,\) by min{/ 7 *(v, \), cy(«, \)} where cy is the capacity function 
of J\f(j'). The modified flow is a preflow for the network M^'). Moreover, d stays a valid 
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Figure 12: A network in which more than one breakpoint is valid. 

distance function as no new residual edges are introduced. Then, we apply the preflow 
algorithm to get a maximum flow /y for the network Af (7'). Since no flow is pushed back 
from the sink, edges from MU2 to \ are kept saturated. Since (7*,k(7*)) and (7', 1(7')) 
are on the same line segment, the minimal cut for /y can also be chosen as (X, X'), which 
indicates that /y saturates all edges to MU\. This implies that 7' is valid for for 
[71,7*). □ 

In Example 15.51 only one breakpoint is valid. In the following example we show that it 
is in general possible that more than one breakpoint is valid. 

Example 5.7. Consider the network depicted in Figure [T2l By a similar analysis as Ex- 
ample [53J we can compute that there are three breakpoints i, 1 and 2. Assuming that 
MU\ = {o\} and MU2 = {02}, we show that all 7 € [53 1] are valid. We send j amount 
of flow along the path /, o\, 02, V, and \ — j amount of flow along the path f,oi,V2,\. If 
7 € [5, 1], then < \ — ? < % implying that the flow on edge fa, \ ) satisfies the capacity 
constraints. Obviously this flow is feasible, and all 7 € [|, 1] are valid for Af(j). 

As we would expect now, it is sufficient to consider only the breakpoints of A/" (7): 

Lemma 5.8. There exists a valid 7 for M{^) iff one of the breakpoints of A/" (7) is valid. 

Proof. If there exists a valid 7 for A/"(7), Lemma [5.61 guarantees that one of the breakpoints 
of ^(7) is valid. The other direction is trivial. □ 

For a given breakpoint, we need to solve one feasible flow problem to check whether it 
is valid. In the network M{^f) the capacities of the edges leading to the sink are increasing 
functions of a real- valued parameter 7. If we reverse A/" (7), we get a parametric network 
that satisfies the conditions in [18]: The capacities emanating from / are non-decreasing 
functions of 7. So we can apply the breakpoint algorithm [18] to obtain all of the breakpoints 
ofA%). 

5.1.3. The Algorithm for DTMCs. Let M be a DTMC and let SimRel„,(A/() denote the 
weak simulation algorithm, which is obtained by replacing line 11.61 of SimRel s (A4) in 
Algorithm [1] by: if s\ t6r S2- The condition si ^.r S2 is checked in Ws(A4, si, S2, R), shown 
as Algorithm 

The first part of the algorithm is the preprocessing part. Line 18.11 tests for the case 
that si could perform only stutter steps with respect to the current relation R. If line 18.51 
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Ws(M,s 1 ,s 2 ,R) 

8.1: if post(si) C i? _1 (s 2 ) then 
8.2: return true 

8.3: if post(s 2 ) C i?(si) then 
8.4: C/i <- {si G posi(si) | si G" i? _1 (s 2 )} 

8.5: return (Vui G Z7l. 3s € post(reach(s2) D R(si)). s G i?(iti)) 

8.6: Compute all of the breakpoints b\ < b% < . . . < bj of J\f (7) 

8.7: return (3i G {1, . . . , j}. 6j is valid for M{b%)) 

Algorithm 8: Algorithm to check whether s\ ^r s 2 . 

is reached, s\ has at least one visible step, and all successors of s 2 can simulate si up 
to the current relation R. In this case we need to check the reachability Condition [3] 
of Definition 12.211 which is performed in line 18.51 Recall that reach(s) denotes the set 
of states that are reachable from s with positive probability. If the algorithm does not 
terminate in the preprocessing part, the breakpoints of the network Affa) are computed. 
Then, corresponding to Lemma 15.81 we check whether one of the breakpoints is valid. We 
show the correctness of the algorithm Ws: 

Lemma 5.9. The algorithm Ws(A4, si, s 2 , R) returns true iff s± ^.r s 2 . 

Proof. We first show the only if direction. Assume that Ws(.M, si, s 2 , R) returns true. We 
consider three possible cases: 

• The algorithm returns true at line 18.21 For this case we have that post(si) C i?~ 1 (s 2 ). 
We choose U\ = 0, V\ = post(si), J7 2 = posi(s 2 ) and V 2 = to fulfill the conditions in 
Definition 12.211 Hence, s\ ^r s 2 . 

• The algorithm returns true at line 18.51 If the algorithm reaches line 18.51 the following 
conditions hold: there exists a state s' x G post(si) such that s[ G" i? _1 (s 2 ) (line 18. 1| ). and 
post{s2) C R{s\) (line I8.3p . Let U\ = {s[ G post(si) \ s[ G" i? _1 (s 2 )}, and define 5{ by: 
<5i(s) = 1 if s G Ui, and otherwise, <5 2 (s) = for all s G S. By construction, to show 
s i ~R s 2 we only need to show the reachability condition. Since the algorithm returns 
true at line !8.51 it holds that for each u\ G Ui, there exists s G post(reach(s2)r\R(s\)) such 
that s G R(ui). This is exactly the reachability condition required by weak simulation 
up to R, thus si ^r s 2 . 

• The algorithm returns true at line 18.71 Thus, there exists breakpoint bi which is valid for 
M(bi). Then, there exists a state s[ G post(s\) such that s[ G" i? (s 2 ), and s' 2 G posi(s 2 ) 
such that s' 2 G" -R(si). By Lemma 15.21 we have that si ^r s 2 . 

Now we show the if direction. Assume that Ws returns false. It is sufficient to show that 
•51 $2- We consider two cases: 

• The algorithm returns false at line l8.51 This implies that there exists a state s[ G post(si) 
such that s[ G" i? _1 (s 2 ) (line 18. ip . and post(s 2 ) C R(s\) (line l8.3p . All states s[ G post(s\) 
with Si G" i? _1 (s 2 ) must be put into U\. However, since the algorithm returns false at 
line 18.51 it holds that there exists a state u\ G U\, such that there does not exist s G 
post{reach(s2)T\R{si)) with s G R(u\). Thus the reachability condition of Definition ^. 211 
is violated which implies that si ^6r s 2 . 

• The algorithm returns false at line 18.71 Then, there exists a state s[ G post(si) such that 
Si G" i? _1 (s 2 ), and s 2 G post(s 2 ) such that s 2 G" -R(si). Moreover, for all breakpoints b of 
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Figure 13: A simple DTMC for illustrating that not all maximum flows are valid. 

Af(j), b is not valid for M(b). By Lemma [5THI there does not exist a valid 7 for Af(-y). By 
Lemma I5T21 we have that si ^.r S2- □ 

Now we state the correctness of the algorithm SimRel^ for DTMCs: 

Theorem 5.10. If SimRel„,(7W) terminates, the returned relation R equals 

Proof. The proof follows exactly the lines of the proof of Theorem 14.21 Let iteration k be 
the last iteration of Ws. Then, by Lemma 15.91 for each pair (si,S2) £ Rk, it holds that 
S2 weakly simulates s\ up to R^, so R^ is a weak simulation. On the other hand, one can 
prove by induction that each i?j is coarser than ^. □ 



Complexity. For (si,sz) € R, we have shown that to check whether s\ tSr S2 we could 
first compute the breakpoints of f\f{"f), then solve C(|V|) many maximum flow problems. 
To achieve a better bound, we first prove that applying a binary search method over the 
breakpoints, we only need to consider C(log |V|) breakpoints, and thus solve C(log|V|) 
maximum flow problems. 

Assume that the sets MU% and PVi for i = 1,2 are constructed as before for N(p(). 
Recall that a flow function / of A/"(7) is valid for N{i) iff / saturates all edges (/, u{) with 
U\ € MU\ and all edges (HJj \ ) with U2 £ MU2- If / is also a maximum flow, we say that 
/ is a valid maximum flow of A/*(7). We first reformulate Lemma 15.21 using maximum flow. 

Lemma 5.11. There exists a valid flow f for N^) iff there exists a valid maximum flow 
f m forM(i). 

Proof. Assume there exists a valid flow / for M (7). We let Aff (7) denote the residual 
network. We use the augmenting path algorithm to get a maximum flow /' in the residual 
network Aff (7). Assume that the augmenting path contains no cycles, which is a harmless 
restriction. Let f m = f + f ■ Obviously, f m is a maximum flow for Af("f), and it saturates 
all of the edges saturated by /. Hence, f m is valid for ^(7). The other direction is simple, 
since a valid maximum flow is also a valid flow for ^(7). Q 

We first discuss how to get a valid maximum flow provided that 7 is valid. Observe 
that even if 7 is valid for Af(j), not all maximum flows for Af(^) are necessarily valid. 
Consider the DTMC on the left part of Figure [TBI Assume that the relation R is given by 
{(si,S2), (si,V2), (vi,S2), (ui,U2), (wi,^)} and consider the pair (si,^)- Thus, we have 
that PVi = {v 1 } 1 MU 1 = {ui},PV 2 = {v 2 },MU2 = {u 2 }- The network ^(7) is depicted 
on the right part of Figure [13j The maximum flow / for AA(1) has value 0.5. If / contains 
positive sub-flow along the path f,ui,V2,\, it does not saturate the edge (H2, \). On the 
contrary, the flow along the single path /, u\,U2, \ with value 0.5 would be a valid maximum 
flow. This example gives us the intuition to use the augmenting path through edges between 
MUi and MU2 as much as possible. For this purpose we define a cost function cost from 
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edges in Af(j) to real numbers as follows: cost{ui,U2) = 2 for u\ G MU\ and u 2 G MU2, 
cost(ui,V2) = 1 for u\ G MU\ and vi G PV2, cost{v\,U2) = 1 for v\ G PV\ and U2 G MU2, 
cost(s,s') = otherwise. The costs of edges starting from source, or ending at sink, or 
in PVi x PV2 are 0. The cost of a flow / is defined by cost(f) = Yle&E /(e) cost (e). By 
definition of the cost function, we have the property cost(f) = f(/~,MUi) + f(MU2,\), i.e., 
the cost equals the sum of the amount of flow from / into MU\ and from MU2 into \,. 

Lemma 5.12. Assume that 7 > is valid for Af (7). Let / 7 denote a maximum flow over 
A/"(7) with maximum cost. Then, / 7 is valid for Nip/). 

Proof. By Lemma I5.11| provided 7 is valid for M{y) , there exists a valid maximum flow 
function g for Mi^y). Since g saturates edges to MU\ and from MU2, obviously, cost(g) = 
P(si, MU\) + 7P(s2, MUz)- Assume that / 7 is not valid, which indicates that / 7 does not 
saturate an edge (/, u\) with u\ G MU\ or an edge (II2, \) with ^2 G MC/2- Then, cost(f 1 ) = 
f{;,MUi) + f(MU2,\) < P(si,Mf7i) +7P(s 2 ,MJ7 2 ) = cos£( 5 ). This contradicts the 
assumption that / 7 has maximum cost. □ 

Let Mu{y) denote the subnetwork of M{j) where the set of vertices is restricted to 
MUi, MU2 and {/A}- The following lemma discusses how to construct a maximum flow 
with maximum cost. 

Lemma 5.13. Assume that f* is an arbitrary maximum flow ofNu{"l) an d f * s an arbitrary 
maximum flow in the residual network A//* (7) with the residual edges from MU± back to f 
removed, as well as the residual edges from \ back to MU2 ■ Then f^ = f* + fisa maximum 
flow over Af{^y) with maximum cost. 

Proof. Recall that the cost of / 7 is equal to cost(f-y) = / 7 (/,M[/i) + f-y(MU2, \). Assume 
that cost(f^) is not maximal for the sake of contradiction. Let / be a maximum flow 
such that cost(f-y) < cost(f). Without loss of generality, we assume that / T (/,M[/i) < 
f{f,MU\). It holds that / 7 = /* + / where /* is a maximum flow of Mu{y), and / is a 
maximum flow in the residual network A//* (7) with the residual edges from MU\ back to 
/ removed, as well as the residual edges from \ back to MU2- On the one hand, /* sends 
as much flow as possible along MU% in A/[/(7). Since in the residual network A//* (7) edges 
from MU% back to / are removed, this guarantees that no flow can be sent back to / from 
MU\. On the other hand, / sends as much flow as possible from MU\ to PV2 (and also 
from PV\ to MU2) in A//* (7). Thus, / 7 (/,MC/ 1 ) must be maximal which contradicts the 
assumption / 7 (/, MU X ) < /(/, MUi). □ 

Assume that 7* is not valid. The following lemma determines, provided a valid 7 exists, 
whether it is greater or smaller than 7*: 

Lemma 5.14. Let 7* G [0, 00) , and let f be a maximum flow function with maximum cost 
for M{y*), as described in Lemma \5.1SX Then, 

(1) If f saturates all edges with u\ G MU\ and (u2,\) with U2 G MU2, 7* is valid 
./'"•AV.'i. ' * _ 

(2) Assume that 3u\ G MU\ such that is not saturated by f, and all edges (u2,\) 
with 1*2 G MU2 are saturated by f . Then, 7* is not valid. If there exists a valid 7, 
7 > 7 . 

(3) Assume that all edges (f~,ui) with u\ G MU\ are saturated by f , and 3u2 G MU2 such 
that (II2, \) is not saturated by f . Then, 7* is not valid. If there exists a valid 7, 7 < 7*. 
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(4) Assume that 3ui € MU\ and 3u2 £ MU2 such that (/,tti) and (U2, \) are not saturated 
by f . Then, there does not exist a valid 7. 

Proof. 1 : Follows directly from the definition. 2 : In this case, < P(si,tti) for 

some u\ 6 MU\. To saturate without un-saturating other edges from / to MUi, 

we have to increase the capacities of edges leading to \, thus increase 7*. 3 : Similar to the 
previous case. 4 : Combining 2 and 3. □ 

According to the above lemma, we can use the binary search method over the break- 
points to check whether there exists a valid breakpoint for ^(7). Since there are at most 
breakpoints, we invoke the maximum flow algorithm at most C(log \V\) times where 
|^| is the number of vertices of A/"(7). 

Theorem 5.15. The algorithm SimRel u ,(A / J) runs in time C(m 2 n 3 ) and in space 0(n 2 ). 
If the fanout g is bounded by a constant, the time complexity is 0(n 5 ). 

Proof. First, we consider a pair (si, S2) out of the current relation Look at a single call 
of Ws(A4, si, S2, Ri). By saving the current relation sets R and in a two dimensional 
array, the conditions s € R(s') or s £ R^ 1 (s') can be checked in constant time. Hence, 
line 18.11 takes time \post(s\)\. To construct the set reach(s) for a state s, BFS can be used, 
which has complexity 0(m). The size of the set reach(s) is bounded by n. Therefore, we 
need O{\post{s\) \ n) time at lines I8.3H8.5I 

The algorithm computes all breakpoints of the network M(^y) (with respect to si,S2 
and R) using the breakpoint algorithm (THJ p. 37-42]. Assume the set of vertices of 7V(7) 
is partitioned into subsets V\ and V2 similar to the network described in Section 14.11 The 
number of edges of the network is at most \E\ := |Vi| | V2 1 — 1- Let \V\ := \V\\ + | V2I , 
and, without loss of generality, we assume that | Vi | < | V2 1 . For our bipartite networks, the 

time complexity [18} p. 42] for computing the breakpoints is C(|Vi| \E\ log(^gy- + 2)) which 

can be simplified to 0(|I4| 2 | V2 1 ) . Then, the binary search can be applied over all of the 
breakpoints to check whether at least one breakpoint is valid, for which we need to solve 
at most Oflog \V\) many maximum flow problems. For our network M{^f), the best known 
complexitju of the maximum flow problem is C(|X^| 3 / log \V\) [2]. As indicated in the 
proof of Lemma 14.41 the distance function is bounded by 4|Vi| for our bipartite network. 
Applying this fact in the complexity analysis in |T3], we get the corresponding complex- 
ity for computing maximum flow for bipartite networks 0(| Vi | \V\ / log | V|). Hence, the 
complexity for the 0(log |V|)-invocations of the maximum flow algorithm is bounded by 
0(|Vi| \V\ Z ). As \V\ < 2\V 2 \, the complexity is equal to \V 2 \ ). Summing over all 

(si,S2) over all Ri, we get the overall complexity of SimRel„,(.M): 

k 

^2 ^2 (\post{si)\ +m+ \post(s 1 )\n+ \Vi\ \V 2 \ 2 ) < Aknm 2 (5.1) 

Recall that in the algorithm SimRel U) (A1), the number of iterations k is at most n 2 . 
Hence, the time complexity amounts to 0(m 2 n 3 ). The space complexity is 0(n 2 ) because 
of the representation of R. If the fanout is bounded by a constant g, we have m < gn, and 
get the complexity C(n 5 ). □ 

^For a network G — (V,E) with small \E\, there are more efficient algorithms in with complexity 
0(\V\ 2 y/\E~\), and in [28] with complexity 0(\E\ \V\ + \V\ 2+€ ) where e is an arbitrary constant. In our 
bipartite networks, however, is in the order of \V\ . Hence, these complexities become 0(11^1 )• 
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5.1.4. An Improvement. The algorithm Ws(A / l, si , S2, R) is dominated by the part in which 
all breakpoints (0(n) many) must be computed, and a binary search is applied to the 
breakpoints, with O(logn) many feasible flow problems. In this section we discuss how to 
achieve an improved algorithm if the network N(^) can be partitioned into sub-networks. 

Let 77 denote the sub-relation R n [(post(si) U {si}) X (post(s2) U {52})], which is 
the local fragment of the relation R. Now let Ai, A2, ■ ■ ■ A^ enumerate the classes of the 
equivalence relation (77 U 77 -1 )* generated by 77, where h denotes the number of classes. 
W. 1. o. g., we assume in the following that A^ is the equivalence class containing si and S2, 
i.e., s±,S2 £ A^. The following lemma gives some properties of the sets Ai provided that 
si £r s 2 : 

Lemma 5.16. For (s\,S2) £ R, assume that there exists a state s^ £ post(si) such that 
s[ £ R~ 1 (s2), and s' 2 £ posi(s 2 ) such that s' 2 £" R(s\). Let A\, . . . ,Ah be the sets constructed 
for (si,S2) as above. If si tSr S2, the following hold: 

(1) P(si, Ai) > and P(s 2 , A) > for all i < h, 

(2) ji = Ki/K 2 for alli<h where ji = P(s 1 ,A i )/P(s2,A i ) 

Proof. Since si ^,r S2, we let 6i,Ui,Vi, A (for 7 = 1,2) as described in Definition 12.211 
Because of states s[ and s' 2 , we have K\ > 0, 7<2 > 0. Let post^Sj) = Ai n post(sj) for 
i = 1, . . . , h and j = 1,2. We prove the first part. For i < h, the set Ai does not contain si 
nor S2, but only (some of) their successors, so it is impossible that both P(si, Ai) = and 
P(s2,Ai) = 0. W. 1. o. g., assume that P(si,Aj) > 0. There exists t £ post^si) such that 
P(si,t) > 0. Obviously 5i(t) = 1, thus: < P(s 1} t) = (K 1 A(t,U 2 ))/5 1 (t) = K 1 A(t,U 2 ) 
which implies that 3u2 £ Ai with A(t, U2) > 0. Hence, P(s2,U2) = ^2A(f7i, it 2 ) > 
7\2A(i, U2) > 0. Now we prove the second part. It holds that: 

P(s 1 ,A i ) = 2^ P(ai,ai) = 2^ P(ai,Oi) = }^ s^{a~) 

{ ^Ki- Y bfaM ^K x - Y, A{a i ,A i ) = K l - Y ^ A i) 

where (*) follows from Condition I2bl of Definition 12.211 (!) follows from the equation <5i(aj) = 
1 for all ai £ post^ai) with i < n, and (f) follows from the fact that if a £ post^si), then 
A(a,b) = for b £ U2\post i (s2). In the same way, we get P(s2,Ai) = K2'Yl a gA dj). 
Therefore, 7, = K1/K2 for 1 < i < h. □ 

For the case h > 1, the above lemma allows to check whether si S2 efficiently. For 
this case we replace lines I8.6H8.7I of Ws by the sub-algorithm WsImproved in Algorithm [9j 
The partition A\ , . . . , Ah is constructed in line 19.11 Lines 19.21 19.101 follow directly from 
Lemma 15.161 if 7$ 7^ 7j for some i,j < h, we conclude from Lemma 15.161 that s\ ^r 52- 
Line 19.111 follows from the following lemma, which is the counterpart of Lemma 15.2b 

Lemma 5.17. For (s\,S2) £ R, assume that there exists a state s'i £ post(s\) such that 
s[ £" R~ 1 (s2), and s' 2 £ post(s2) such that s' 2 £" R(si). Assume that h > 1, and assume 
WsImproved(.M, si, S2, R) reaches line \9.11\ Then, s\ T^r S2 iff 71 is valid for N (71) . 

Proof. First, assume that s\ ^r S2- According to Lemma 15.21 there exists a valid 7* for 
A/"(7*). As in the proof of Lemma 15.21 7* = K1/K2 is valid for M(j*). If Ws reaches 
line 19.111 by Lemma 15.161 we have 71 = K1/K2, hence, 71 is valid for AT(7i). The other 
direction follows directly from Lemma 15.21 □ 
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WsImproved(A1, si, s 2 , R) 

Construct the partition Ai, ... , Ah (* Assume that h > 1 

for all i <— 1, 2, ... h — 1 do 

if P(si,A) =P(s 2 ,^i) = 0then 

raise error 
else if P(si, Aj) = or P(s 2 ,Aj) = then 
return false 

else 

P(siA) 

if 7j ^ 7j for some i, j < /i then 

return false 
return (71 is valid for A^(7i)) 

Algorithm 9: Algorithm to check whether s\ Tin s 2 tailored to DTMCs. 



Example 5.18. Consider the DTMC in Figure [8] together with the relation R = {(si,s 2 ), 
(s 1 ,v 2 ), (vi,s 2 ), (ui,u 2 ), (01, o 2 ), (01, v 2 ), (vi,o 3 ), (vi,v 2 ), (o 2 , 0^}. We obtain the relation 
H = R\{(o 2 , 01)}. We get two partitions A\ = {ui,u 2 } and A 2 = {s%, s 2 , vi, v 2 , 01, o 2 , 03}. 
In this case we have h = 2. Recall that MU\ = {ui,o\}, MU 2 = {u 2 , 02,03}, PV\ = {^i}, 
PV 2 = {v 2 }. WehaveP(si,Ai) = \ and P(s 2 , Ax) = |. Hence, 71 =P(s 1 ,A 1 )/P(s 2 ,A 1 ) = 
2. As we have shown in Example 15.41 2 is valid for the network M{2). Hence, s\ t6r s 2 . 

Assume that (si, s 2 ) G R\ such that h > 1 in the first iteration of SimRel^. We consider 
the set A\ and let 71 = P(si, Ai)/P(s2, Ai). If A\ is not split in the next iteration, 71 would 
not change, and hence, we can reuse the network constructed in the last iteration. Assume 
that in the next iteration A\ is split into two sets A\ and A\. There are two possibilities: 

• either P(s 1; Af)/P(s 2 ,A\) = P(s 1 ,Af)/P(s 2 ,A\). This implies that both of them are 
equal to 71. If all A4 are split like A±, we just check whether 71 is valid for M^i). 

• or P(s l5 A%)/P(s 2 , A\) / P(si, Af)/P(s 2 , A\). This case is simple, we conclude s 1 £ R s 2 
because of Lemma 15.161 

This indicates that once in the first iteration 71 is determined for (si, S2), either it does not 
change throughout the iterations, or we conclude that s\ s 2 directly. The above analysis 
can be generalised to the case in which Ai is split into more than two sets. As the network 
Af("fi) is fixed, we can apply an algorithm similar to Smf, which solves the maximum flow 
problems during all subsequent iterations using only one parametric maximum flow, as for 
strong simulation. 

The above analysis implies that if h > 1 for all (s\, s 2 ) in the initial R\, we could even 
establish the time bound 0(m 2 n), the same as for strong simulation. Since in the worst 
case it could be the case that h = 1 for all (si,s 2 ) £ R, the algorithm WsImproved does 
not improve the worst case complexity. 

Since the case that the network cannot be partitioned {h = 1) is the one that requires 
most of our attention, we suggest a heuristic approach that can reduce the number of 
occurrences of this case. We may choose to run some iterations incompletely (as long as 
the last iteration is run completely). If iteration i is incomplete, we first check for each pair 
(si,s 2 ) £ Ri whether the corresponding hi is greater than 1. If not, we skip the test and 
just add (si,s 2 ) to Ri+i- The intuition is that in the next complete iteration i' > i, for 
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each such pair (si, S2) we hope to get > 1 because some other elements of Ri have been 
eliminated from it. We only perform the expensive computation if an incomplete iteration 
does no longer refine the relation. 

5.2. An Algorithm for CTMCs. Let M = (S,K,L) be a CTMC. We now discuss 
how to handle CTMCs. Recall that in Definition 12.231 we have the rate condition [3j: 
Ki~R(si, S) < K2R(s2, S). To determine ^m, we simplify the algorithm for DTMCs. If 
K\ > and K2 = 0, s± ^.r S2 because of the rate condition. Hence, we do not need to 
check the reachability condition, and lines I8.3ll8.5l of the algorithm Ws(M, si, S2, R) can 
be skipped. For states si, S2 and relation R, we use ^(7) to denote the network defined in 
the embedded DTMC emb(Ai). To check the additional rate condition we use the following 
lemma: 

Lemma 5.19. Let s\ R S2- Assume that there exists s^ € post(s\) such that s\ ^ R~ 1 (s2). 
Then, s± S2 in M. iff there exists 7 < R(s2, S*)/R(si, S) such that 7 is valid for Mi^y). 

Proof. Assume first s± S2 in M.. Let 5i, Ui, Vi, Ki, A (for i = 1,2) as described in 
Definition 12.231 Obviously, s^ must be in U\, implying that K\ > 0. Because of the rate 
condition it holds that KiR(sx, S) < l^R^j S), which implies that K2 > 0. It is sufficient 
to show that 7 := K1/K2 is valid for M{y). Exactly as in the proof of Lemma 15.21 (the 
only if direction), we can construct a valid flow / for J\f (7). Thus, 7 is valid for J\f (7) and 
7<R(s 2 ,«S)/R(si,S). 

Now we show the other direction. By assumption, 7 is valid for M{y). We may assume 
that there exists a valid flow function / for M (7). We define <5j, Vi, Ui, Ki, A (for i = 1,2) as 
in the proof (the if direction) of Lemma 15.21 Recall that s^ must be in U\, implying that 
/(/,s'i) > 0. Thus, there must be a node s in M(y) with f(s,\) > 0, which implies that 
s € U2. Thus we have K2 > 0. Using the proof (the if direction) of Lemma 15.21 it holds 
that s\ ^Sr S2 in emb(A4), moreover, it holds that 7 = K\/K~2. By assumption it holds that 
7 < R(s2, 5')/R('Si, S) which is exactly the rate condition. □ 

To check the rate condition for the case h > 1, we replace line 19. Ill of the algorithm 
WsImproved by: 

return (71 < 7* A 71 is valid for N{y\)) 
where 7* = R(s2, 5)/R(si, S) can be computed directly. In case h = 1, we replace line 18.71 
of Ws by: 

return (3i G {1, . . . , j}. hi < 7* A bi is valid for N(bi)) 

to check the rate condition. Or, equivalently, we can check whether the minimal valid 
breakpoint y m is smaller than or equal to 7*. The binary search algorithm introduced for 
DTMCs can also be modified slightly to find the minimal valid breakpoint. The idea is 
that, if we find a valid breakpoint, we first save it, and then continue the binary search on 
the left side. If another breakpoint is valid, we save the smaller one. As the check for the 
reachability condition disappears for CTMCs, we get better bound for sparse CTMCs: 

Theorem 5.20. If the fanout g of M is bounded by a constant, the time complexity for 
CTMC is 0(n 4 ). 

Proof. In the proof of Theorem 1 5 . 1 5 1 we have shown that Ws has complexity C?(|Vi| | V2 1 2 ). 
As we do not need to check the reachability condition, the overall complexity of the al- 
gorithm SimRel^M) (see Inequality El]) is E Sl es E S2 eS £i=i (\post( Sl )\ + |Vi| \V 2 \ 2 ) 
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which is bounded by 2kgm 2 . Since k is bounded by n 2 , the time complexity is bounded by 
Agm 2 n 2 . If g is a constant, we have m < gn, hence, the time complexity is 4g 3 n 4 G C(n 4 ). □ 

6. Conclusion and Future Work 

In this paper we have proposed efficient algorithms deciding simulation on Markov models 
with complexity 0{m 2 n). For sparse models where the fanout is bounded by a constant, 
we achieve for strong simulation the complexity 0{n 2 ) and for weak simulation 0(n 4 ) on 
CTMCs and 0(n 5 ) for DTMCs. We extended the algorithms for computing simulation 
preorders to handle PAs with the complexity 0(m 2 n) for strong simulation. For strong 
probabilistic simulation, we have shown that the preorder can be determined by solving 
LP problems. We also considered their continuous-time analogon, CPAs, and arrived at an 
algorithm with same complexities as for PAs. 
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